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Abstract 

This thesis develops a method to analyze the maneuvering forces on surfaced and 
underwater vehicles with complex propulsors. The analysis method is developed for 
general propellers yet has unique applicability to model highly contracting stern flows 
associated with integrated propulsors. Integrated propulsors exhibit strong coupling 
of the various blade-rows and duct, if present, to the vehicle stern. The method 
developed herein provides a robust means to analyze propulsor-induced maneuvering 
forces including those arising from wake-adapted, multi-stage, ducted propulsors. 

The heart of the maneuvering force prediction is a three-dimensional, unsteady 
lifting-surface method developed as the first part of this thesis. The new method is 
designated PUF-14 for Propeller Unsteady Forces. The lifting-surface method uses 
many advanced techniques. One significant advance is the use of a wake-adapted 
lattice to model the flow through the propulsor. In related research, a 2-D Kutta 
condition has been augmented using Lagrangian interpolation to dramatically reduce 
the required computational time to model a 2-D gust. 

The second thrust of this thesis couples the unsteady lifting-surface method with a 
three-dimensional, time-average Reynolds-Averaged Navier-Stokes flow solver. Rotat- 
ing a propeller through a spatially-varying flow field causes temporally-varying forces 
on the propeller. From the converged-coupled solution, the maneuvering and blade- 
rate forces can be estimated. This thesis explores the relationship of time-varying 
and time-average forces in the flow solver and potential-flow domains. Similarly, it 
explores the relationship of the effective inflow in the two domains. Finally, this the- 
sis details the synergistic means to correctly couple the potential-flow method to a 
viscous solver. 

Verification and validation of the method have been done on a variety of geometries 
and vehicles. Preliminary results show good correlation with experiment. The results 
strongly suggest this maneuvering force prediction method has great potential for the 
modern propulsor designer. 

Thesis Supervisor: Justin E. Kerwin 
Title: Professor of Naval Architecture 
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Chapter 1 



Introduction 

1.1 Overview 

In order to study complex propulsors, a method of solving the resulting time-dependent 
non-linear boundary value problem is required. Methods of predicting the propeller 
blade-row performance in steady and unsteady flow have been in existence for many 
years. Comprehensive reviews of the steady and unsteady lifting-surface theory are 
given in Kerwin [27] and Schwanecke [34], respectively. 

The present work builds on advances in steady lifting-surface theory, as in PBD- 
14 [28]. and advances in unsteady lifting-surface theory, such as PUF-2 [31]. By 
blending the best features of the two methodologies, a robust methodology is devel- 
oped to calculate the time- varying forces on highly complex propulsor geometries. 
Then, the new lifting-surface methodology is mated with a three-dimensional, steady 
Reynolds- Averaged Navier-Stokes (RANS) solver. The RANS solver captures the 
viscous nature of the propulsor characteristics and its interaction on the body. The 
resulting coupled methodology, documented in [42], is able to analyze the steady and 
unsteady blade-row forces on highly complex propulsors. 

1.2 Need for a New Methodology 

Figures 1-1 and 1-2 illustrate examples of highly complex propulsors. The hydrody- 
namic complexity arises due to strong interaction of the rotor with the full afterbody. 



1C 



The duct, if present, tends to inhibit flow separation as the body contracts at the 
vehicle stern. The propulsor components are strongly coupled and must be designed 
and analyzed while taking into account their full interaction. This type of propulsor 
- where strong interaction occurs between components and the body - is known as an 
integrated propulsor. The behavior of an integrated propulsor during a maneuver is 
difficult to predict with cylindrical-streamtube design methods. The current research 
should provide a relatively fast and robust analysis method for integrated propulsors. 

Figure 1-3 shows the typical body dynamics during a maneuver. During a given 
phase of the turn, the body sees an inflow that varies slowly compared to the propeller 
rotation. This inflow can be represented with a time-average RAXS solution and 
the propulsor interaction modeled using the new methodology. Figures 1-4 and 1- 
5 illustrate the body-centered coordinate system used when describing maneuvering 
forces and moments. 




Figure 1-1; Notional submerged-body integrated propulsor (Developed at MIT). 



Methodologies to calculate unsteady forces on highly complex propulsors, such as 
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Figure 1-2: Notional surface ship underwater hull form with an integrated propulsor concept (De- 
veloped at NSWC under the electric drive program.) 

figure 1-1, exist today only as rudimentary approximations. The new methodology 
provides a method to assess the unsteady forces arising in the wake-deficient regions 
behind control surfaces, etc. and those arising due to once-per-revolution spatial 
variations as the vehicle maneuvers. In the end, the coupled methodology should be 
invaluable to the modern propulsor designer. 

1.3 Functional Description of the Thesis 

The areas of concentration of this thesis work are divided into two functional groups. 
First, the vortex-lattice lifting-surface methodology is developed. An original lifting- 
surface program, PUF-14, has been written to support the new methodology. This 
first functional group develops a method to perform an analysis of the blade-row 
in a specified inflow. Areas of advancement include: wake-adaptive modeling in an 
unsteady method, incorporation of hub and duct images in an unsteady method, 
and Glauret spaced panels which led to a very promising “de-singularized" Kutta 
condition. 

Second, the lifting-surface method is coupled with a RANS code. Our ONR spon- 
sor has provided a three-dimensional, steady RANS code. The RANS code is used 
to obtain the inflow velocity field, which the lifting-surface code uses to calculate the 
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Figure 1-3: Typical maneuvering transient. Courtesy of SNAME. 

unsteady forces. Unsteady forces are generated due to rotating the propeller in a 
spatially-varying inflow. Time-averaged, but spatially-varying body forces are intro- 
duced into a three-dimensional volume to represent propulsor stages in the RANS 
flow field. The entire RANS flow field responds to the blade-row presence. In turn, 
the RANS flow field is used again for the lifting-surface analysis of the blade-row. By 
alternately updating the lifting-surface and the RANS solutions, the blade-row forces 
and RANS flow field converge to the appropriate solution. Areas of advancement 
include: determination of the proper velocities to obtain the correct volumetric effec- 
tive inflow, determination of the proper velocities to obtain the correct time-average 
forces and detailed study of the relationship of time-average velocities and forces in 
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Figure 1-4: Typical steady inflow during maneuvering transient is at an angle /3 with respect to the 
body centerline. 
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Figure 1-5: Typical steady inflow during depth transient is at an angle a with respect to the body 
centerline. 

both the RANS and the potential-flow domains. 

The new treatment of unsteady force calculations should greatly improve propulsor 
prediction capabilities. In practice, past and current-day efforts treat unsteady forces 
with simplified assumptions such as cylindrical propulsor geometry and extremely 
limited body /propulsor interaction, if any at all. Competing with the new coupled 
approach are fully three-dimensional, unsteady RANS formulations, which are being 
developed at other institutions. While unsteady RANS offers great potential, the 
computing burden is enormous and not yet practical in modern applications. The new 
treatment developed in this thesis is believed to be practical in both computational 
load and in representing the physical hydrodynamic characteristics of today’s complex 
propulsors. 
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Chapter 2 



Formulation of the Lifting-Surface 
Problem 

2.1 Fundamental Assumptions 

The propeller is assumed to be a set of thin blades arranged symmetrically about 
a common axis. No restriction is made on the blade shape. The propeller rotates 
with constant angular velocity in an unbounded fluid. Hub and duct effects are 
represented with blade images. The onset flow is permitted to be a function of all 
spatial coordinates. 

The fluid is assumed to be incompressible and the flow field irrotational 1 except 
on the blade and in the trailing vortex wake sheets. The blade boundary layer and 
shed vortex wake thickness is assumed to be thin so that the fluid rotation due to the 
propeller is confined to a thin layer. The fluid is treated as inviscid except for some 
empirical corrections associated with the propeller blade drag. 

2.2 Boundary Value Problem 

In the fluid domain, the velocity, V, must satisfy 

V V = 0 (2.1) 

l To be precise, the prescribed inflow field may be rotational, but vortical interaction with the potential flow field 
induced by the blade singularities is not treated explicitly in the blade solution, but rather, is accounted for by 
interaction with a RANS code. 
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( 2 . 2 ) 



or. with the assumption of the previous section, 

V 2 d> = 0 

where 4> is the velocity potential defined by \ = V0. In a propeller-fixed reference 
frame, 

V -72 = 0 (2.3) 

on the camber surface. At large distances from the propeller 

V -> U (2.4) 

where l is the specified onset flow. 

Additionally, Kelvin's theorem for the conservation of circulation is explicitly in- 
voked on the blade and in the wake sheet. The Kutta condition is implicitly satisfied 
on the trailing edge of each blade. Further, since a wake sheet is unable to support 
a pressure jump across it. the vorticity in the wake must align with the velocities in 
the wake to remain force free. 

2,3 Singularity Distribution 

The solution of the boundary value problem is approached by distributing singularit ies 
on the mean camber surface of the propeller and on the vortex wake. Thus, the field 
equation, V 2 <£ = 0, is immediately satisfied. 

A distribution of sources on the camber surface is used to generated the jump in 
normal velocity. The source strengths are obtained by stripwise application of thin 
wing theory. The jump in tangential velocity across the camber surface is provided 
by a distribution of vorticity on the camber surface and in the wake sheet. The 
vorticity strengths are obtained by imposition of the boundary conditions stated 
above. This leads to a singular integral equation over the blade and wake surface. By 
the numerical scheme described later, the boundary problem is reduced into a system 
of linear algebraic equations. 
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Chapter 3 



Mathematical Modeling of the 
Propeller System 

3.1 Blade Geometry 

Propeller blades are traditionally defined by tabular data consisting of radial distribu- 
tions of geometrical quantities such as pitch, rake and skew, and by chord wise distri- 
butions of camber and thickness defined along circular cylinders. However, problems 
arise due to interpolation and smoothing during each step in the design and man- 
ufacturing process. This problem is intensified as the design process becomes more 
complex, where geometry must be passed among a large number of hydrodynamic 
and structural analysis codes. 

To uniquely define the blade geometry, B-spline surfaces are used. A B-spline 
surface can be interrogated at an arbitrarily fine mesh of points to define the blade. 
Therefore, designers and manufacturers each can use the same B-spline to define the 
geometry to the required accuracy of their respective calculations [35]. 

The propeller blade surface is described by means of a control polygon net of B- 
spline vertices. The number of control polygon vertices required to define the camber 
surface of a propeller blade can be quite small, with a 7 x 7 grid being satisfactory 
in many applications. The B-spline surface is easy to manipulate during the design 
process, such as when designing with PBD-14 [39], because each B-spline vertex has 
a local effect as shown in figure 3-1. 
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Control Net and Surface Gnd 



Figure 3-1: The effect on blade shape of moving one B-spline control vertex. 

Pie- and post-processors enable the use of the designer's description of the blade 
surface, i.e pitch, rake, skew. The designer's descriptions are necessary when commu- 
nicating with older codes/tools. The tools used for these conversions are described in 
the Propeller Blade Design (PBD-14) Manual [39]. 

The blade coordinate system used in this work closely follows that described in 
[29] and [14]. A Cartesian coordinate system is fixed to the propeller with the x- 
axis coaxial with the propeller hub and pointing in the streamwise direction. The 
y-axis may be oriented at any convenient angle to the propeller. The corresponding 
cylindrical coordinates are defined with r 2 = y 1 + z 1 and the blade rotation angle, 6. 
is measured from the y-axis in a right-handed sense. Figure 3-2 illustrates the blade 
geometry and the coordinate system. 

3.2 Discretization of Blade Singularity Distribution 

The method of singularity distribution is one of the most powerful techniques for 
the solution of the fluid flow problem. The boundary value problem formulated in 
chapter 2 can be divided into two distinct groups: the source singularities and the 
vortex singularities. 
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Figure H-*2: Pictorial of the blade coordinate system. Note that a right-handed propeller rotates in 
a direction of negative angle. 

The continuous distribution of sources and vortices is replaced by an array of 
M x iX concentrated straight-line elements of constant strength. The end points 
of the elements are located on the camber surface. The exact arrangement of the 
points has been the subject of much experimentation. The chosen spacing on the 
camber surface follows the method elected by Greeley and Kenvin [14]. Figure 3-3 
illustrates the straight-line representation of the continuous singularity distributions 
for a M x A r = 11 x 5 blade grid. 

3.2.1 Source Distribution 

rhe source singularities are distributed to represent the jump in normal velocity across 
the camber surface. At the outset, we assume the source strength distribution is 
independent of time, and that its spatial distribution may be derived from a stripwise 
application of thin-wing theory at each radius. These assumptions are consistent with 
Kenvin and Lee [29]. 
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Figure 3-3: This figure shows the blade grid and wake grid points connected by a mesh. The 
blade control points are also shown. The boundary value problem solves for the blade circulation 
distribution which simultaneously nulls the normal velocity at every control point on every blade 
surface. 



Strictly speaking, the source strengths are dependent on the local inflow and will 
therefore vary in time. However, by using the mean inflow to set the source strengths, 
the Variations with time from the mean value will be small. Additionally, the blade 
thickness contribution to the boundary value problem is secondary considering both 
the mean and fluctuating blade loading. In recognition of the secondary influence 
of blade thickness, it is a reasonable assumption that the source strengths are time 
invariant. Thus, the source strength and their contribution to the boundary value 
problem are solved only once. 

3.2.2 Vortex Distribution 

The jump in t angential velocity across the camber surface is provided by a dist ribut ion 
of vorticity on the camber surface and in the wake sheet. As shown in figure 3-3, the 
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vortex strength is a vector lying on the surface and may be resolved into components 
along two arbitrarily assigned directions on the surface. The vortex distribution on 
the blade is chosen to be resolved into "spanwise” and “chordwise” components while 
the corresponding components in the wake are termed ‘‘shed” and "trailing" vorticity. 
In general, the shed vorticity is a result of temporal variations of the spanwise vorticity 
strength on the blade. 

An implicit Kutta condition is established by arranging the blade grid and the 
control points such that a control point lies on the blade trailing edge. Thus, when 
the condition of V’ -ft = 0 is satisfied exactly the blade trailing edge, the Kutta 
condition is simultaneously satisfied. The implicit Kutta condition has been shown 
to work in Greeley and Kerwin [14] and Breslin. et al. [4]. 

3.3 Geometry of Transition and Ultimate Wakes 

The geometry of the trailing vortex wake has an important influence on the accuracy 
of the calculation of induced velocities on the blade. The wake alignment procedure 
employed by Greeley and Kerwin [14] convects the trailing vortices along axisymmet- 
ric surfaces constructed from a user-supplied tip contraction angle and ultimate wake 
radius. The present scheme follows that used in PBD-14 and allows the wake to con- 
form to the actual inflow velocity field. Thus, the wake conforms to the body shape 
and interior duct surface, if present, and follows the circumferential-mean stream- 
tubes. 

The wake is divided into two regions: 

1. Transition Wake. The transition wake region, wherein all contraction and 
roll up of the trailing vortex sheet is assumed to occur, is built upon the actual 
inflow velocity field. The transition wake is an extension of the blade vortex- 
lattice grid. The axial extent of the transition wake is specified by the user. 
The shed vorticity in the transition wake is assumed to decay in strength as it 
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is washed downstream. 



2. Ultimate Wake. The ultimate wake comprises the region from the end of the 
transition wake to a point infinitely far downstream. The trailing vortex lines 
from each blade are merged into an infinite-bladed helical vortex whose radii 
match that of the downstream end of the trailing vortex lines. The total cir- 
culation of the helical vortex is determined by the average of the shed vorticitv 
at each individual trailing vortex. Efficient closed-form expressions for the ve- 
locity induced by infinite-bladed helical vortices were developed by Hough and 
Ordwav [15]. Leibman [32] demonstrated that the infinite-bladed approximation 
introduced negligible errors provided that the transition wake length was at least 
one propeller radius. 

The transition wake grid is implicitly specified by the grid density on the blade 
and by the time stepping variables. The trailing vortices in the transition wake are 
extensions of the chordwise vortices on the blade. The discretized representation of 
the vortex sheet then consists of .1/ concentrated trailing vortex lines whose trailing 
edge coordinates match the corresponding values of the chordwise vortices on the 
blade. 

The transition wake region contains a set of A j x M shed vortex segments, where 
Aj is the number of shed vortex segments in the transition wake. The segments 
occupy the region from the blade trailing edge to the start of the ultimate wake 
The total inflow which convects the transition wake elements can vary with angular 
position as seen in figure 3-4. However, as a simplifying assumption, the wake is built 
upon the circumferential-mean inflow. The wake geometry has no dependence on the 
blade angular position; thus, the wake geometry remains unchanged relative to the 
blade. Thereby, a single wake geometry can be used to represent all transition wakes 
at all angular positions. While this is artificial, it is deemed necessary considering 
the scope of the project and is expected to add a negligible error in the calculations. 



28 




Figure 3-4: This figure shows a center-body representation and the blade grid overlayed with the 
input axial velocity. 

Tracking of individual wakes trajectories could be added in a future revision. 

The angular increment, SO , and the time increment, St , are related through the 
equation 

SO = uoSt (3*1) 

where w is the rotation rate of the propeller. 

Given the non-dimensional assumptions of Vs = 1 and Rbiadetip = 1 , the convection 
parameter which governs the time interval between each time step is 

* = t < 3 - 2 > 

where Js is the advance coefficient and is the number of time steps per revolution. 
Thus, the convection of the shed vorticity is related to all components of the inflow 
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velocity and to the time increment implicitly specified by the input parameters. This 
improves upon the linear convection used in Kerwin and Lee [29] and Keenan [24]. 
[23]. ■ Additionally, by design, there exists an explicit correspondence of the shed 
vortex geometry to the time step interval, the n th shed vortex at time t will have the 
same strength as the (n — l) </l vortex at time (t — 1). 

While the shed vorticity is convecting downstream, it is hypothesized that the 
shed vorticity becomes highly disorganized within a short distance behind the trailing 
edge [29]. In order to simulate this dissipation of shed vorticity, a decay factor 7 (P) 
is applied to the shed vorticity strength given by 

7(P) = 2P 3 - 3P 2 + 1 0 < P < 1 (3.3) 

where P is the fractional distance along the vortex sheet from the blade trailing edge 
to the ultimate wake [23] . In this way, the shed vortex strength is made to approach 
the average value for that particular M set of vortices. Conveniently, the ultimate 
wake approximation uses this average value to estimate the strength of the ultimate 
wake. Thereby, a smooth transition of vortex strength is maintained from the blade 
through the transition wake and into the ultimate wake. 

3.4 Representation of Propeller and Wake Vorticity 

The propeller is assumed to have several blades which are identical in shape and 
evenly spaced around a common axis. Some earlier analysis techniques solved the 
boundary value problem by solving the problem on only one blade. 1 However, for 
this solution method, the boundary value problem includes all blades. Thus, every 
control point on every blade is simultaneously solved. 

The effects of the vortex elements are taken into account by conceptualizing the 
vorticity as closed loops which automatically satisfy Kelvin’s theorem. Figure 3-5 

1 Specifically, PBD-14 is an steady, axisymmetric solver which solves one blade BVP. The influence of all blades are 

correctly accounted for through the influence function. Conversely, PUF-2 used a coarse spacing on all blades except 

the key blade and iteratively approached a solution at each time step. 
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shows the conceptual idealization. The influence of a blade vortex element on the 
velocity at a control point is captured with a Loop Influence Function. The induced 
velocities due to the LIF are unknown for the current time step. The LIF is placed 
on the left-hand-side of equation 3.4. 

The influence of a wake vortex element is captured using a Wake Influence Func- 
tion , WIF . The wake influence depends on the time history of blade loading, thus the 
induced velocities on a control point are known quantities. Consequently, the wake 
vorticity-induced velocities are placed on the right hand side of equation 3.4. 

The LIFs and WIFs are developed as the velocity induced on a control point due 
to a unit strength vortex loop. Lee [31] and Keenan [24] show the construction of the 
simultaneous equations for closed loops. The present work follows Lee and Keenan 
in the derivation of loops. 

Lhilike Lee and Keenan, the present work solves the boundary value problem on 
all blades simultaneously. The formulation of the simultaneous equations, while con- 
ceptually similar, yields many more equations. An additional difference is that the 
present formulation uses cosine spacing across the chord of the blade with an implicit 
Kutta condition at the trailing edge. 

As expected, the LIF matrix exhibits the structure for the propeller system. Equa- 
tions 3.4 and 3.5 show the structure that exists in the simultaneous equations. These 
equations show the influence functions after they have been converted to represent 
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Figure 3-5: The boundary value problem is solved by conceptualizing the vorticitv as closed loops 
on the blade and in the wake. The ultimate wake is a “closed” loop extending to infinity. 



the normal component of velocity due to a unit-strength vortex loop. 
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(3.4) 

Accounting for symmetry such as Blade-on-Blade which is the same for all blades 
and Blade-on-Prior-Blade which is the same for all Blade-on-Prior-Blade relations. 
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then the structure is more evident. Equation 3.5 assumes a five-bladed propeller and 
replaces LIF with .4. LIF'j 2 with B. etc., to use the Blade-on-- • • relations which 
makes the structure much more clear. 

Vortex Elements 
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3.5 Stepwise Solution in the Time Domain 

With the knowledge of how the propeller system is to be modeled, the sequence of 
program execution must be considered. Figure 3-6 shows a functional block diagram 
of how the program execution should proceed. Note that the simultaneous equations 
must be solved at each time step. All information that is not time dependent is 
calculated prior to entering into the time domain calculations. 
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PUF14 Flow Chart 




Figure 3-6: Functional block diagram of the PUF-14 methodology 
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Chapter 4 

Determination of Blade Forces 



The force and moment acting on the propeller blade can be obtained by integrating 
the pressure jump over the blade camber surface, 

F = I A/md.4 (4.1) 

M = J Ap(r x n )<LA (4.2) 

where S]> is the pressure jump across the camber surface, n is the normal vector 
on the camber surface, and r is the position vector from the origin to the point of 
integration. 

The forces acting on a propeller blade can be divided into four portions: the 
pressure forces acting normal to the blade surface, the viscous forces acting tangential 
to the blade surface, the leading edge suction force and the force proportional to the 
time rate of change of potential, which follows from the unsteady term in Bernoulli's 
equation. The four components of blade force are discussed below. 

4.1 Pressure Forces 

The force acting normal to the blade surface, commonly referred to as the pressure 
force, is obtained from Joukowski's law. At the control point in a given panel there 
is a vortex density 7 corresponding to the local pressure jump across the blade. The 
Joukowski force on the panel is 

dFj — p dA (V x 7 ), (4.3) 



35 



where V is the mean velocity acting at the control point of the panel. dA is the area 
of the panel and p is the fluid density. The vortex density 7 is related to the jump in 
tangential velocity, 2 c/v, across the blade by 

7 = n x 2 dv, (4.4) 

n being the unit normal to the blade surface. This can be used to rewrite the 
Joukowski force as 

dFj = pdA ((V • 2dv) n - (V • n)2dv) . (4.5) 

The jump in tangential velocity, 2 e/v, is obtained by taking the gradient of the 
potential jump across the vortex sheet 

2 r/v = Vs//. (4.6) 

Vs is the gradient operator on the blade surface and p is the potential jump across 
the vortex sheet. In a vortex-lattice representation, the potential jump at the control 
point is equal to the vortex loop strength around the panel, so the jump velocity is 
obtained directly by differentiating the vortex loop strengths representing the loaded 
blade surface. 

4.2 Viscous Forces 

Viscous forces are computed using an empirical sectional drag coefficient. This coef- 
ficient, Cdv, is supplied by the user. An element of viscous force acting on one panel 
of the blade surface is then computed as 

dF v = l -pV\V\C dv dA . (4.7) 

The total drag force is then obtained by summing the elemental forces over all panels. 

4.3 Leading Edge Suction Forces 

A blade which is not acting at ideal angle of attack (shock-free entry) has an additional 
force acting on it — the leading edge suction — which is a consequence of representing 
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the blade as a thin wing. The suction force on a segment of the leading edge is 
computed as 

dF L = -j~pC*b , (4.8) 

where b is a unit vector, orthogonal to the leading edge, lying on the camber surface 
of the blade and C s is the suction coefficient. This coefficient is obtained from the 
limit 

C s = lim \/s ( 7 ( 5 ) • t), (4.9) 

s— >0 

where t is a unit vector tangent to the leading edge, and s is the curvilinear distance 
along the camber surface from the leading edge. The total suction force is obtained 
by summing the elemental dF^s over the length of the leading edge [25]. 

4.4 Unsteady Velocity Potential 

As the blade passes through a region where the inflow changes, the blade experiences 
an added-mass force. This force is proportional to the time derivative of the veloc- 
ity potential and follows from the unsteady term, equation 4.10, from Bernoulli's 
equation. 

dF A = pn-^(4> + - 4)~)dA (4.10) 

First, recall that the source strengths are assumed to the independent of time, so that 
only the potential due to the vortices is required. Since the jump in potential across 
the camber surface in the continuous case is 

4>+(s)-4>-(s)= [ s 7(CMC (4.11) 

J L.E. 

we can replace 4.10 by the final form 

dF A = pn-j- [ i(()d(dA (4.12) 

at J l.e. 

The time derivative in 4.12 may be obtained by numerical differentiation of the dis- 
crete vortex strengths obtained at several discrete time steps. Since the calculations 



37 



of the forces can be made after the results of the propeller have been time-domain 
solution, the future is already known: thus, permitting the use of a very accurate 
central five-point differentiation formula. 
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Chapter 5 



Validation using a Specified 
Effective Inflow 

5.1 Convergence 

There are many parameters that must converge to yield the final converged answer. 
The first parameter studied is the convergence of the solution for various blade grid 
resolutions. Propeller 4577 grid resolution is varied from a coarse resolution of S x 8 
lattice to a much finer lattice of 36 x 36. The chosen inflow is described in section 5.3. 
The results are plotted in figure 5-1. The figure shows good convergence of I\ j and 
I\q for a relatively coarse grid. In practice, the desired grid is one which gives the 
desired accuracy with the least computational effort. A grid of around 12 x 12 provides 
a reasonable compromise. 

The next convergence study centers on unsteady analysis. Propeller 4118 is used to 
test the results of PUF-14 for convergence in an unsteady analysis. The chosen wake is 
the 4118 wake, published in [29], with slight modifications at the hub and tip to ensure 
the wake extends to the blade extremes. The wake is essentially described as one plane 
of flow data which varies circumferentially and with radius. This particular 4118 wake 
was experimentally made by superimposing wakes that have strong harmonic content 
in the third and fourth harmonic. 

The solved circulation for the entire time domain solution is shown in figure 5-2. 
When the '‘starting vortex” is first generated, there is a strong effect on the solved 
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Figure 5-1: Convergence of t he forces with increasing blade grid resolution. 

circulation. As the “starting vortex" washes downstream, its strong local effects are 
mitigated which allows the solved circulation to approach its converged value for each 
particular angular position. 

The next significant parameter to test for convergence is the size of the time step. 
This parameter, called the number of time steps per revolution or A#, controls the 
size of the angular rotation per time step which, in turn, controls the downstream 
convection distances of each vortex loop. Figure 5-3 shows the error associated with 
various time step sizes. The plotted parameter is the circulation at the 0.7 r/R of the 
blade. The circulation is strongly influenced by the harmonics in the inflow velocity 
field. The converged solution for circulation using 210 time steps per revolution is 
taken as the reference circulation. Thus, the plotted values represent the percent 
error associated with a given Nq. 

Note that in figure 5-3 the average magnitude of error does not improve signifi- 
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Specified Tolerance = 0.0001 
Max actual (GAM at 0.7)=0.0495 
2.00 rev Max tol. = 0.00932 
2.45 rev Max tol. = 0.0000891 




Figure 5-2: The solved circulation is shown for the entire time domain solution calculated by Pl'F-14. 
Note that only the converged solution correctly represents the boundary value problem. 

cantlv above 60 X$. Thus, a value of 60 A,;, is a good compromise between speed 
and accuracy. Figure 5-4 shows a direct comparison of two values of the Xg. 60 and 
210. These two cases show that the solved blade circulation tracks closely for the two 
values of Ng . 

In viewing figure 5-3, there is obviously a periodicity to the error associated with 
the convergence of the time step size. The maximum error tends to occur where the 
solved circulation has the greatest slope. Therefore, a distinction in the error associ- 
ated with amplitude and phase may be seen through a comparison of the harmonic 
analysis. 

Figures 5-5 and 5-6 show the decomposition of the convergence error associated 
with the force calculations into amplitude and phase. The amplitude is fairly well 
converged. However, the phase is not so settled and, in fact, becomes cpiite poor 
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Figure 5-3: Convergence error associated with the angular size of the time step, X# 




Figure 5-4: Comparison of two angular sizes of the time step, N$ 



in the higher harmonics. The errors in phase are somewhat expected due to the 
arbitrary choices in converting the continuous boundary value problem into a discrete 
problem of simultaneous equations. The arbitrary choices become less important as 



the discrete time step becomes smaller to approach the continuous case. Therefore, 
when accurate phase information is need at the higher harmonics, smaller time steps 
are needed. 




Figure 5-5: Comparison of harmonic phases of various time step sizes, A'^ 




Figure 5-6: Comparison of harmonic amplitudes of various time step sizes, 
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5.2 Validation Against Other Numerical Methods 



The stand-alone mode of PUF-14 was validated in an effective wake by comparing to 
the results of PUF-2.1 and PBD-14.2.S in the identical wakes, respectively. Since these 
two codes. PUF-2.1 and PBD-14.2.S, have been in existence for many years, they have 
been validated both at the initial code validation phase and, later, validated against 
real-life designs. Thus, the PUF-14 computer program can be validated against these 
two pre-existing programs for the specialized cases discussed below. 

5.3 Steady Analysis 

PBD-14. 2. 8 has been validated in such references as [11]. [2], [1]. Comparison against 
PBD-14 will be used to validate PUF-14 in the stead}' mode of analysis. The chosen 
test case uses the propeller 4577 from Huang [IS]. The effective inflow is approximated 
by using the total inflow from a converged axisymmetric RANS flow field for the 
Huang body one. While the effective inflow may not be correct for this propeller, it 
provides a rough estimate of an effective flow field and allows a direct comparison of 
the results of the two codes. 

5.3.1 Solved Blade Circulation 

The two cases shown in figure 5-7 illustrate the agreement between PUF-14 and PBD- 
14. 2. S. Additionally, cases run using various combinations of parameters, including 
stator blade-rows, all show nearly exact agreement between PBD-14. 2. 8 and PUF-14 
for steady axisymmetric test cases. 

5.3.2 Forces 

The two cases shown in figure 5-7 have a steady value of thrust and torque. These 
values are shown in tables 5.1 and 5.2. The results are exact, indicating that PUF-14 
is correctly solving the special case of steady, axisymmetric propeller analysis. As 
stated above, additional cases all show similar agreement. 



44 



Propeller 4577 in a notional 




Figure 5-7: Comparison of the Non-Dimensional Circulation versus Radius 
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PBD- 14.2.8 
PUF-14 


0.0838 

0.0838 


0.1920 

0.1920 



Table 5.1: Steady axisymmetric force comparison for the case with thickness, hub, duct. G blades 
with a 10x10 grid 





At 


Kq 


PBD-14.2.8 

PUF-ll 


0.0567 

0.0567 


0.1219 

0.1219 



Table 5.2: Steady axisymmetric force comparison for the case with no hub, duct nor thickness, 3 
blades with a 12x12 grid 



5.4 Unsteady Analysis 

PUF-2.1 has been validated by many authors [31], [29], [26]. Thus, comparisons 
with PUF-14 will suffices for validation of correctly solving the specialized problem 
of circumferential- and radial-varying inflow but no variation along the propeller axis 
and no hub nor duct. The chosen test case is propeller 4118 as discussed earlier in 
this chapter. 
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The most significant difference between the codes leading to the different solved 
circulation is the lattice arrangement on the blades and in the wake. Figure o-S shows 
the key blade and the transition wake for the two codes. Recall that PUF-14 uses 
identical blade lattices and identical transition wakes for all blades. PUF-2.1 uses a 
much coarser lattice on blades other than the key blade. 




Figure 5-b: The leftmost figure shows PUF-2.1 key blade and transition wake lattice arrangement. 
The rightmost figure shows PUF-14 key blade and transition wake lattice arrangement. Differences 
in lattices lead to differences in solved blade circulation shown later. 

More specifically. PUF-2.1 uses constant spacing in the chord wise direction while 
PUF-14 uses cosine spacing. Therefore, PUF-14 resolves more of the leading and 
trailing edge gradients. Another lattice difference is in the wake models: PUF-2.1 
only propagates about a quarter of a propeller diameter aft while PUF-14 propagates a 
user-specified distance aft. Another difference exists in the specification of the general 
shape for the transition wake and in the ultimate wakes. Still another difference 
exists in the specification of the blade geometry and lattice. That is, PUF-2.1 uses 
cylindrical geometry parameters to built its lattice while PUF-14 uses a B-spline 
geometry description overlayed with a vortex lattice built on the circumferential- 
mean streamtubes. 
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5.4.1 Solved Blade Circulation 



Figure 5-9 shows the agreement of the blade circulation at the 0.7 r/R position. Some 
differences of circulation are expected due to the inherent differences in the solution 
methodology between the two codes. 




Figure 5-9: Comparison of the Non-Dimensional Circulation at 0.7 r/R 

5.4.2 Forces 

File cases shown in figures 5-10 and 5-11 compare the x, y , and z forces and mo- 
ments, respectively, on a blade as it rotates. Note that the xyz coordinate system 
rotates with the blade. As discussed above, exact agreement of the solved circulation 
and forces is not expected. In the figures, there is evidence of a small error in the 
magnitude of forces/moments and a small error in the phase of the forces/moments. 
However, there is significant agreement between the two solution methods which leads 
to the conclusion that PUF-14 is correctly solving the specialized, unsteady propeller 
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boundary value problem of cylindrical flow. 



5.5 Validation Against Experiment 

5.5.1 Experimental Setup 

Experiments are described in reference [21] in which the propeller 4679 is used to 
obtain unsteady pressure distributions in inclined flow. The propeller was designed 
to model a controllable-pitch (CP) propeller used for high-powered ships. A picture 
of propeller 4679 is shown in figure 5-12. 

The experiment was conducted on the DTXSRDC Carriage Y. The propeller was 
driven from downstream while the carriage supporting the propeller advanced into 
the flow. The comparisons chosen for validation of PUF-14 use the test in which the 
propeller shaft was inclined at 7.5 degrees downward from the direction of the carriage 
advance so that the propeller operated in inclined flow. Two advance coefficients have 
been used for comparison of PUF-14 and this experiment. 

Xo flow field measurements were given so the geometric flow field is used to es- 
timate the propeller inflow. The inclined flow is represented by once-per-revolution 
variation in the radial and tangential flow components. As a further refinement of the 
flow field, the prescribed geometrical inflow was modified with a crude adjustment 
for the presence of the propeller shaft. An axial potential flow approximation was 
included by modeling the propeller shaft as a Rankine ovoid centered on the propeller 
plane. A similar correction was applied for the cross-flow directions using the poten- 
tial flow around a two-dimensional cylinder. Convection velocities were used for wake 
alignment as computed from PSF-2 [14]. 

5.5.2 Comparisons with Experiment 

Figures 5-13 and 5-14 show the mean pressure distribution on one blade of DTMB 
4679 at r/R = 0.7 for Js of 1.078 (design Js) and 0.719, respectively. The figures 
compare PUF-14, PUF-10.3 and the experimental results. PUF-10 is an unsteady, 
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low-order perturbation-potential-based panel method incorporating an unsteady it- 
erative pressure Kutta condition at the blade trailing edge. The trailing vortex wake 
has a frozen geometry based upon user supplied inflow inputs and convection veloc- 
ities. PUF-10 was given the same specified inputs as PUF-14. PUF-10 provides a 
check of the trends expected for the pressure distribution. The experimental results 
provide point values expected from the calculations. PUF-10. 3 and PUF-14 follow 
the same basis shape while nearly matching the experimental data points. 

Figures 5-15 through 5-18 show a comparison of the first harmonic amplitudes 
and phases. The figures compare the experimental results with the calculational 
results of PUF-14. The PUF-10 results are not in agreement with either PUF-14 
or the experiment so those results are not shown here. 1 A definite agreement is 
established between PUF-14 calculations and the experiment. The two agree for on- 
and off-design advance coefficients in both the first harmonic amplitude and the first 
harmonic phase. 

The reasons for the differences that do exist are unclear. One reason is likely due 
to the source strength representation of blade thickness. In PUF-14. the line sources 
used to represent the blade thickness are solved at the onset of the solution using the 
circumferential-mean inflow. Therefore, the source strengths are invariant in time and 
space. This assumption means that the once-per-revolution variation in tangential 
velocity on the blade, which linearly determines the source strengths, is not captured. 
Thus, the net impact is that velocities variations at higher harmonics are not totally 
correct. That is, the velocity does not include the perturbation velocities which would 
be due to the source strength variation as the blade rotates. This simplification 
does not significantly affect the force calculations since forces are derived from a 
pressure difference. However, the simplification is likely to affect, to a slight degree, 
the calculated pressure distribution mean value as well as its harmonic variations. 

1 It is likely that the unsteady case was incorrectly setup for the PUF-10 calculations. As such, this statement is 
not meant to judge the capability of the PUF-10 panel method. 
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Comparisons at r/R = 0.5 and r/R = 0.9 radii were also made for the mean, first 
harmonic amplitude and first harmonic phase. The comparisons show similar agree- 
ment between PUF-14 and the experiment. The comparisons with the experiment 
confirm that PUF-14 is correctly solving the boundary value problem for this case of 
inclined effective inflow using modestly skewed propeller blades. 

5.5.3 Propeller AY and I\q 

Finally, in the effort to compare with this experiment, many methods were used 
to validate the various aspects of the calculations. Thus, as a by-product of this 
experimental comparison, the propeller forces and moments were calculated by several 
different methods. Figure 5-19 presents the calculations for the mean values of AY 
and I\ Q at the design conditions of the propeller. The comparisons of the mean value 
of AY and I\q suggest that PUF-14 is correctly solving the global blade-row forces 
and moments, but may be over predicting the mean values by a few percent. 
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: igure 5-10: Comparison of the PUF-14 and PUF-2.1 force analysis 
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Figure 5-11: Comparison of the PUF-14 and PUF-2.1 moment analysis 
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Figure 5-12: Propeller 4679 blade shape for the unsteady calculations. 



52 




Figure 5-13: Mean pressure distribution on one blade of DTMB 4679 at r/R — 0.7 for the design J$ 
of 1.078. 




Figure 5-14: Mean pressure distribution on one blade of DTMB 4679 at r/R — 0.7 for Js of 0.719. 
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lire 5-15: First harmonic amplitude pressure distribution on one blade of DTMB 4679 at r/R 
for the design Js- 




Figure 5-16: First harmonic amplitude pressure distribution on one blade of DTMB 4679 at r/R 
0.7 for Js of 0.719. 
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Figure 5-17: First harmonic phase pressure distribution on one blade of DTMB 4679 at r/R =0.7 
for (he design J 5 . 




Figure 5-18: First harmonic phase pressure distribution on one blade of DTMB 4679 at r/R = 0.7 
for J s of 0.719. 
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Figure 5-19: Propeller 4679 mean A't and I\q for the design J$ = 1.078 and an off-design Js = 0.719. 
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Chapter 6 



Coupling with Three-Dimensional, 
Steady Viscous Flow Solver 

6.1 Overview of Coupling 

The concept of coupling a potential flow propeller unsteady forces (PUF) analysis 
method with a viscous flow solver is relatively simple. The difficulty of designing 
coupled methods comes in resolving the details of the processes. In practice, the 
coupling of lifting-surface theory with viscous flow solvers has proven useful in ax- 
isymmetric cases. Some early efforts of coupling potential flow analysis methods 
with viscous solvers were by Science Applications International Corporation [43] and 
Kerwin, et.ai [2S] in 1993 and 1994, respectively. 

Prior to the implementation of coupling methodologies, perhaps the most trou- 
blesome assumption had been that the propeller operates in potential flow. In fact, 
propellers almost always operate in the highly rotational shear flows of the ship's 
boundary layer and wake. The propeller does not operate in the flow field measured 
in nominal wake surveys and one designed to that flow would often be deficient in 
actual use. This circumstance leads to the so-called effective wake , or effective in- 
flow, problem. The need for an effective inflow description arises from the flow at 
the propeller's location being different with and without the propeller present. This 
difference is due to coupling of the propeller’s induced velocity field to the vorticity 
in the incoming flow. A redistribution of the vorticity in the inflow, and hence the 
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shear profile of the flow at the propeller, results. 

Potential flow theory provides a powerful basis for representing the propeller's 
own flow field but contains no mechanism to treat the effective inflow problem. On 
the other hand, viscous flow methods naturally capture the vortical phenomena of 
the inflow to the propeller but offer a poor framework for representing the propeller 
itself. The coupled method described herein couples the two formulations, viscous 
flow and potential flow, to rationally address the effective inflow problem in a fully 
three-dimensional flow field. 

To handle the effective inflow problem requires a flow solver that can explicitly 
model the transport of vorticity. There are a number of methods that do this. An 
additional requirement, though, is that the solver be able to capture separation. 
Integrated-stern concepts involve hull forms that actually depend on the action of 
the propulsor to inhibit separation as discussed in Dai et.al. [9], Chen tt.nl. [5]. 
Warren [41], Thus, a viscous flow solver, specifically, a Reynolds-Averaged Xavier- 
Stokes (RAXS) solver is used. 

The method of using a lifting-surface with a viscous solver has the significant 
advantage of easily supporting multiple blade-row analysis. To treat such a problem in 
potential flow alone leads to numerical difficulties as wake-sheet singularities approach 
control points on the downstream blade row. In the coupled method, all the vorticity 
is dealt with by the RANS code so that there are no singular structures and the 
velocity field is smooth. In addition, the potential-flow propeller analysis only has to 
deal with one blade-row at a time. 

The first part of the problem treats the viscous flow around a body while the second 
part treats the inviscid problem of the flow around the blade-row. The presence 
of the blade-row in the viscous solution is represented by a suitable distribution 
of body forces obtained from the blade-row analysis. In return, the viscous flow 
solver provides a distribution of total velocity which serves as the input to the blade- 
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row analysis. However, the precise relationship between body forces and blade-row 
inflow velocity is not obvious. In particular, this portion of the thesis explores the 
extraction of the effective inflow from the total inflow in fully three-dimensional flow 
and the relationship of the forces to the total and effective inflow. Additionally, this 
thesis explores the relationship of the flows during three-dimensional coupling when 
including hub/duct images, multiple blade-rows and a highly contracting flow field. 

Recently, coupling a PUF code with a RANS code was done at the Naval Surface 
Warfare Center, Carderock Division as a proof of concept for the current work [19]. 
Their coupling successfully showed the concept of three-dimensional coupling for an 
open-water single blade-row case. Although, their coupling used a three-dimensional 
RANS code, it only extracted a two-dimensional inflow for computing the effective 
wake. 1 

6.2 Relationship of Time- Average and Unsteady 

Imagine a fully-appended underwater body being towed forward through the water. 
From the vantage point of the body, the gross flow field does not change with the 
advancing of time. There will be some turbulent eddies which seemly fluctuate ran- 
domly with time, but the magnitudes of the fluctuations are small compared to the 
gross flow field velocity. The flow field can be averaged in time to remove the turbu- 
lent fluctuations. Then, from the vantage point of the body, the flow field is constant 
in time - or in the language of this thesis - the flow field is time- averaged. 

The appendages leave a momentum defect associated with their wakes. Towards 
the body stern, there are regions near the freestream velocity, regions of higher velocity 
and regions of lower velocity. By using process of time averaging, the flow field retains 
the flow variations around the body’s features. 

J The NSWC coupling used DTNS3D and PUF-2. PUF-2 assumes that the velocity just upstream of the propeller 
can be used everywhere over the axial extent of the propeller. Its formulation ignores axial variations in the inflow 
velocity. Thus, the NSWC coupling used a two-dimensional effective inflow disk at a predetermined axial location 
just upstream of the propeller and ignored axial variation and radial contraction of the inflow. 
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Now imagine the underwater body has a propeller attached to it’s stern which 
pushes the body through the water. At a fixed point in space near the propeller, the 
passage of the propeller blades produces a time- varying velocity. As shown some time 
ago by Breslin [4], this fluctuating velocity field is due to a combination of the motion 
of the blades relative to the observer, and the time-varying loading on the blades 
themselves. This gives rise to velocity fluctuations at all harmonic orders of shaft 
rate - not just blade rate. In this way, rotating a propeller through a time-average 
flow field can produce propeller loading which varies in time - or in the language of 
this thesis - produces unsteady propeller loading. 

6.3 Lifting-Surface/RANS Coupling Details 

The lifting-surface method solves the time-varying boundary value problem. These 
solutions give blade loading as a function of time and correspondingly as a function 
of blade angular position. The blade loading is translated into stationary body forces 
in the RAXS domain. The body force influences the RANS solution, which in turn, 
influences the flow field near of the propeller. Then, the flow field is extracted from 
RANS and used as input to the lifting-surface solution. The whole cycle is repeated 
until convergence. Figure 6-1 shows the coupled solution methodology. Figure 6-2 
shows the some of the details of an iteration. The coupling suite of utility programs 
are labeled in these figures. Their names are: cs_setup - setup program, cs.unns - 
velocity extraction program, and cs.chipr - force interpolation program. 

6.3.1 Time- Average Induced Velocity 

Some means must be formulated to get the effective inflow from the flow solver. The 
goal is to extract the time-average effective inflow from the time-average total inflow 
as determined by the viscous flow solver. The relationship of effective inflow to the 
total inflow is given in equation 6.1. The axisymmetric effective inflow is discussed 
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Figure 6-1: Coupled solution methodology flowchart 
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Figure 6-2: Coupled solution iteration details 
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in numerous papers including Kenvin et. al. [28]. 
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The time- varying loading on a rotating propeller will induce a time- varying velocity 
at a fixed point in space near the propeller. Some means must be made to average 
the velocity fluctuations due to a rotating propeller. Figure 6-3 shows the calculated 
induced velocities for one fixed point in space due to the rotation of a three-bladed 
propeller. The plotted data describes the time-varying induced velocity produced by 
all blades during one propeller revolution. The figure also shows a constant line at the 
value of the time-average velocity resulting from averaging the time-varying velocity 
produced by the blade-row. This time-averaging is necessary at all field points in the 
swept volume of the propeller. 2 

The current methodology uses the fully three-dimensional swept volume of the 
propeller to couple with the viscous solver. As expected, then, the time-average 
induced velocity must be calculated at every control point at every blade position. 
The entire swept volume is described at field points which have been especially chosen 
to coincide with the control points whenever the blade passes the field point. This 
method removes the singular effect of a vortex approaching too near a field point. 
Additionally, by ensuring the solution time step, A*, is evenly divisible by the number 
of blades, a field point is guaranteed to be co-located with a control point at every 
blade time-step position. 

After obtaining the time-average induced velocity, then the effective inflow can be 
obtained by the relationship in equation 6.1 and seen visually in figure 6-4. Figure 6- 
4 shows some interesting features. The upstream flow has a velocity deficient over 
about SO degrees of the inflow disk. When the propeller passes the region with the 
slower velocity, the blade loading increases. Higher blade loading results in higher 

2 As a computational note, the calculation of the time-average induced velocity is extremely computationally 
intensive. The difficulty lies with the need to capture the large velocities produced by the blade vortex sheet as a 
field point is brought very close to the sheet. This need must be balanced with the need to minimize the required 
computational time. This calculation could be improved in the future. 
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Velocity at FP due to all three blades 





Blade Position (degrees) 

Figure 6-3: This specific example shows the induced velocity at one field point (FP) near mid-span 
and mid-chord of the three-bladed propeller. 

induced velocity as seen by the lighter colors in the upper-right plot. The velocity 
deficient remains in the effective inflow. Note that downstream of the propeller, the 
velocity deficient has causes additional contraction of the propeller tip streamtube 
due the higher blade loading in the slower velocity region. The contraction is seen 
as darker colors in the total inflow plot. The darker colors show the relatively slower 
freestream being pulled in to the plotted flow domain. 

6.3.2 Time-Average Forces 

It is desired to have the viscous flow solver solve for the time-average velocity that 
represents the flow in the presence of a propeller. The forces calculated by the pro- 
peller analysis must somehow be placed into the flow domain. One method introduces 
the time-average forces into the flow domain as body forces . 3 In this way, when fluid 

3 Body forces are analogous to gravity forces. Gravity is a body force that is directed downward and is generally 
constant in strength. A propeller body force is placed in the flow solver with direction and strength according to the 
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Total Inflow - Induced Inflow = Effective inflow 




Figure 6-4: Given the total velocity (from RANS), the time-average induced velocity is subtracted 
away leaving the effective inflow in which the boundary value problem is solved. 



is in the swept volume of the blade-row, the fluid is accelerated according to the forces 
produced by the blade-row. Thus, the fluid feels the presence of the blade-row and 
the viscous solver produces a solution that represents the time-average total inflow to 
the blade-row. In this work, based in part on the conclusions of Stern et. al. [37. 36]. 
blade-row forces are placed in the viscous flow solver as bod)' forces in the swept 
volume of the blade-row. 

Time- averaging of the body forces must be accomplished such that the time- 
average body force is consistent with the flow-solver time-average velocity. Kerwin 
et. al. [28] discusses the relationship of propeller forces to body forces for the circum- 

blade-row solution. 
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ferential mean flow field as well as the relationship of the total circumferential inflow 
to the effective inflow. 

Kerwin et. al. [28] shows that the forces that must be applied in the flow solver to 
produce the correct effective inflow arise for the circumferential mean inflow and the 
blade vorticity. The circumferential mean inflow is distinct for the local inflow. The 
local inflow is the inflow exactly at the blade mean-camber surface. The circumferen- 
tial mean inflow is the inflow that is averaged at one field point during one propeller 
revolution. 

The extension to three-dimensional unsteady analysis is as follows. To place the 
correct forces into the flow solver, the forces must be calculated using the blade 
vorticity and the time-average velocity. That is. to keep the time-average velocities the 
same in the flow solver and in the blade BVP, an equivalent force must be transmitted 
to the flow solver. The equivalent force will maintain equal circulations in the flow 
solver and the blade BVP. Using a force calculated with the blade local velocities 
would be incorrect and lead to erroneous results when coupled with the flow solver. 
See section 7.1 for more discussion on the coupling method, and reference [28] for 
more discussion on the forces required to obtain the correct axisymmetric solution. 

6.3.3 Specific Details of a Rotor 

With each revolution, the rotor blade sweeps a constant volume of fluid. The revo- 
lution is discretized by the number of time steps per revolution, No . such that blade 
forces are calculated at Ng locations. Discretizing in this manner yields the blade 
forces as a function of angular position. The discrete approximation of forces is made 
continuous in the RANS domain by interpolating the forces in the region between 
discrete blade N$ positions. Figure 6-5 shows one meridional plane which contains 
rotor forces interpolated onto the RANS cells. Figure 6-6 show's the forces in the 
swept volume of the rotor. 

Bv interpolating forces between discrete blade positions, the RANS domain is pro- 
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Figure 6-5: This figure shows one meridional-plane containing blade forces interpolated to the viscous 
solver's cells. 



vided a continuous distribution of forces in the swept volume of the rotor outline. 
Additionally, the forces in the swept volume are time-average forces. Local effects 
such as secondary flows and tip vortex flows cannot be captured due to the contin- 
uous nature of the forces. While the mean and unsteady forces are captured, the 
contribution to these forces due to secondary flows is not captured. 



6.3.4 Specific Details of a Non-rotating Blade Row (Stator) 



The stator blades could be gridded with RANS cells just as the body, duct and 
appendages are gridded. This would leave only the rotor to be modeled with the 
body forces in the current methodology. However, the burden of gridding the stator 
details could be removed by modeling the stator with the body forces. Given that 
the rotor is already modeled with the body forces, it seems appropriate to model the 
stator in the same manner. In the work presented herein, the stator is not gridded 
with the RANS grid. 

The representation of a stator in this methodology is similar to the rotor rep- 
resentation with some unique differences. The stator is represented in RANS with 
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Figure 6-6: This figure shows a notional inflow with the corresponding time-average forces in the 
swept volume of the rotor. (Color reproduction in figure B-l.) 



body forces, which are only placed where the mean-camber blade surface overlaps the 
RANS cells. 

The correct forces to represent a stator in RANS are those forces generated with 
time-average induced velocities (TAV) just the same as for a rotor. However, for a 
stator, the TAV are, in fact, identical to the local blade induced velocities (LBV). 
One unique difference between the stator and rotor representation is that the stator 
forces are not distributed in a continuous swept volume like rotor forces. The stator 
forces are local to the blade surface, which does not rotate. Thus, the forces are only 
located in the RANS cells only where the stator surface overlaps the RANS cells. By 
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representing the stator in this way, the rotor analysis will capture some of the stator 
blade-rate harmonic forces due to the stator interaction. 
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Chapter 7 



Preliminary Methodology 
Validation 

7.1 Propeller 4119 

7.1.1 Setup 

The DTMB Propeller 4119 is used here to perform a preliminary validation of the 
coupling procedure. The case is close enough to an ideal case that some aspect can 
be verified by inspecting the results. The DTMB Propeller 4119 is a three-bladed 
propeller of relatively simple geometry. The propeller was designed by Denny for 
uniform inflow as a double thickness version of Propeller 4118 [10]. Experimental data 
taken by Jessup includes measured pressure distributions, boundary-layer profiles, 
momentum based drag coefficients and downstream wake surveys for this propeller 
at the design advance coefficient [22]. Using the new coupling methodology, the 
measured circulation and the experimental results for I\j and I\q are compared to 
the coupled results 

For the water tunnel experiment, the diameter-based Reynolds number (Rep) was 
766,400 at the design advance coefficient of 0.833. The chordwise Reynolds numbers 
varied from 0.8 to 2.0 million from hub to tip. Experimental results were obtained 
with the blade both tripped and smooth. The tripping was done with roughness at 
the leading edge. 
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7.1.2 Calculations 



The propeller is placed in a cylindrical flow domain within the RANS domain. The 
chosen RANS grid is 89 axially by 81 radially by 25 azimuthally. The grid is rather 
coarse but has been proven adequate for these results. Figure 7-1 shows one merid- 
ional plane of the grid and the outline formed by revolving the grid about the x-axis. 
The swept volume of propeller 4119 is shown. The grid extends 4 propeller radii 
upstream. C radii downstream and 12 radii radially. The specified upstream inflow is 
axisymmetric with U = 1 and V = W = 0. 




Figure 7-1: Straight-shaft RANS grid, grid outline, and the swept volume of propeller 4119 

It should be noted that the convergence of the numerics is first checked using a 
stand-alone PUF evaluation in an assumed effective inflow. The blade lattice con- 
verges using a 15x15 lattice to within one percent of the converged value of I\j. The 
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AY converges when using variable Ng within about a very coarse 10 steps per revo- 
lution. However, Ng of 10 is unacceptable when performing coupled analysis. In a 
coupled. evaluation using PUF and a RANS code, the time-average velocities play an 
integral part in ensuring convergence to the correct solution. At present, the calcula- 
tion of the time-average velocities uses Ng angular positions. As Ng, becomes smaller 
the accuracy of the calculation for the time-average velocity becomes much worse. 
Therefore, the conclusion is that Ng should be a minimum of about 60. 

Figure 7-2 shows the numerical dynamics during the coupled analysis. The figure 
shows three different sets of coupling solutions. The set marked with an open circle 
is PBD-14 coupled with the axisvmmetric version of DTNS. The set marked with a 
U is PUF-14 coupled with the three-dimensional serial-UNCLE. The set marked with 
D is PUF-14 tricked into coupling with the axisymmetric version of DTNS. Note 
this is an axisymmetric test case. By using three-dimensional analysis tools, their 
preservation of the axisymmetric properties is check and the actual results can be 
directly compared with the PBD/DTNS coupling solution that has been verified by 
others. There appears to be some error associated with the three-dimensional coupling 
that can be either attributed to the coupling suite of utility codes or attributed to the 
three-dimensional RANS code. For this special case, the three-dimensional coupling 
should closely match the PBD/DTNS coupling. The cause of this three-dimensional 
error is not obvious. 

During a coupled analysis, the numerical dynamics followed the expected trends. 
At iteration 0 (PUF analysis in the nominal inflow of U = 1 and V = W = 0), the AY 
is expected to be too low. The reason is that the nominal inflow does not include the 
induced velocities in the propeller wake. This causes the wake pitch to be too tightly 
coiled, which in turn, causes too much induced velocity at the propeller region. Large 
induced velocity unloads the propeller and causes AY to be too low. 

At the next iteration (the second PUF analysis), the AY is too high. The reason 
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Propeller 4119 




Figure 7-2: Coupled analysis convergence history for propeller 4119 



is as follows. Iteration 0 body forces are low which causes the total velocity from 
RANS to be smaller that it should be. Additionally, Iteration 0 induced velocity, 
which is subtracted from the total, is too large due to the tightly coiled wake. Both 
the total too low and the prior induced too high cause the current effective inflow to 
be too low. This combination results in AY too large. The dynamics will eventually 
converge to the correct solution, with the effective inflow equal to the nominal inflow 
and with the propeller wake aligned with the total flow field. 

One of the most useful aspects of this test case is in comparing the effective inflow 
to the expected effective inflow. The flow has an upstream velocity of U = 1 with 
V = W = 0, i.e. circumferential-mean inflow. No vorticity exist in the flow at the 
upstream boundary. The grid was not intended to be fine enough to accurately cap- 
ture the boundary layer. However, the boundary layer, which does build on the shaft, 
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is confined to a small percentage of the blade span. Additionally, the experimental 
setup uses a shaft with a nose cone upstream of the propeller which changed the 
propeller inflow velocity profiles from the profile that is being calculated here with 
RANS. This effect is neglected in this study since reasonable results can be obtained 
with less computational time using a straight-shaft grid. 

For this special case (of using a straight shaft) away from the shaft boundary layer 
effects, the nominal flow 1 must be equal to the effective inflow. This is true since 
there is no pre-existing vorticity in the inflow to interact with the blade induction. 

A confirmation of the three-dimensional error can be seen in the effective inflow. 
After convergence of the coupled lifting-surface/RANS solution, the effective inflow 
to the blade boundary value problem must be U = 1 with V = IT = 0. Figure 7-3 
shows the effective inflow after the convergence of the lifting-surface/RANS solution. 
Propeller drag and thickness have been set to zero. This check must hold true for any 
propeller, with any number of blades and at any advance coefficient. Many variations 
of the check were tried and each check resulted in similar effective inflow. 

7.1.3 Possible Source of Error 

For completeness, this section is included. Many avenues were tried to uncover why 
the effective inflow does not match the expected inflow. Sources of error include a 
myriad of possibilities. However, one source has stood out. The two RANS codes, 
DTNSax verses UNCLE, provide different velocities when given the identical body 
force distribution. Note that the effective inflow was low by 1-2% when PUF was 
coupled with either the DTNS3D or UNCLE. Leading one to conclude the if the error 
is, in fact, in the RANS code, then the error is in the body force routines (which 
are common) or in a fundamental assumption of the RANS formulation or in some 
combination. 2 

x That is, the flow that would be present in the absence of the propeller. 

2 This section is not drawing the conclusion that the error is in the RANS solvers, but simply stating that the 
difference in the 3-D solvers and the axisymmetric solver is possibly the source of error. More detailed study of this 
difference in velocity is necessary before more definitive conclusions can be drawn. 
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Figure 7-3: Away from hub and tip effects, the axial effective inflow for propeller 4119 on a straight 
shaft should be equal to 1.0. 



The 3-D RANS solution estimates the axial velocity to be a few percent lower 
than the axisymmetric RANS. With a total velocity a few percent lower, the effective 
velocity will be lower. Of course, the non-linear nature of coupling a propeller with 
a RANS code could mitigate or enhance this effect. But, the trend it consistent with 
the straight-shaft effective inflow being too low and consistent with the resulting 
I\'t being too high. Figure 7-4 shows the RANS results for the identical body force 
distribution. The body forces used for this comparison used only axial body forces so 
that any doubts about coordinate transformations between axi- and 3-D coordinate 
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systems are removed. 



Radial and Tangential Velocity 




Figure 7-4: Velocity comparison at the propeller plane (x/R = 0) between DTNSax and UNCLE for 
identical axial-only body force distributions. 



7.1.4 Circulation of Propeller 4119 

Using Laser Doppler Velocimetry (LDV) the flow field was measured at two planes 
downstream of the propeller. The circulation in the flow field was computed by 
Jessup by integrating the tangential velocity at the upstream plane. The experimental 
results are compared to the lifting-surface calculations in figure 7-5. Note that the 
experimental data shown has had the circulation due to the tangential component of 
the boundary-layer wake subtracted. This was done by Jessup to allow comparison 
with potential flow circulations. 

The experimental results are not expected to exactly match the calculated results. 
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Foremost, the RANS domain does not include the nose cone leading the shaft. Sec- 
ondly, the differences in the circulation near the tip are due to the wake contracting 
slightly and a tip vortex beginning to form at the measurement plane. Near the hub, 
the boundary-layer calculations may be under-predicting the flow separation that 
occurs there. 

Figure 7-5 does confirm the that the circulations calculated in the lifting-surface 
problem are very near to the circulations placed in the RANS domain. Using theoret- 
ical extrapolations discussed in section 6.3.2, this figure supports the notion that the 
correct forces are being transferred into the RANS domain. However, this figure also 
indicates that the coupled lifting-surface/RANS solution is resulting in higher blade 
loading than expected, which can be reasoned as consistent with an effective inflow 
that is too low. 

7.1.5 Blade-row Force Prediction 

An additional check with the propeller 4119 is a check of the mean AY and I\q. The 
experimental results at design J=0.833 are AY = 0.145 and I\q = 0.0280. Including 
the effects of thickness and drag, the calculated results are AY = 0.160 and I\q — 
0.0293. Recall the RANS grid does not accurately model the nose cone leading 
the shaft. Some error associated with the mean AY and I\q is expected due to 
the assumptions of the flow domain geometry and some error is associated with the 
incorrect effective inflow. 

7.1.6 Potential Error in the Coupling Mechanics 

The coupling methodology is illustrated in figures 6-1 and 6-2. The weights block 
shows the pre-calculation for the interpolation. The weights dictate the fraction of 
forces that go into specific RANS cells. Pre-calculating weights saves computational 
time. However, if the weights are not updated and the nominal weights are used 
throughout the coupled analysis, then the effective inflow will be in error by about 
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Propeller 4119 




Figure 7-5: Circulation distribution for Propeller 4119 from experiment, PUF-14 lifting-surface so- 
lution and from the RANS solution. 

2% for this straight-shaft case. This 2% error in effective inflow would be added to 
the error previously discussed. For this case, the 2% error in effective inflow will lead 
to about a 3-4% error in I\j and I\q . To remove this possible source of error, the 
weights are updated periodically during a coupled analysis. 

7.2 MIT Pre-swirl 

7.2.1 MIT Pre-swirl Experimental Setup 

In an attempt to validate the representation of the stator when using the coupled 
methodology, the MIT Pre-swirl experiment was used as a validation case. Forces 
on the stator and rotor combination were obtained experimentally in the MIT water 
tunnel. The stator was designed by Bowling [3] to provide swirled inflow to David 
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Taylor Model Basin model propeller 4497. The stator geometry and rotor geometry 
are described in [3]. The experimental geometry, as represented in the coupled analysis 
method, is shown in figure 7-6. The experiment measured mean rotor thrust and 
torque. 




Figure 7-6: Coupled analysis method representation of MIT Pre-swirl geometry. Flow would be 
generally from right to left in this figure. 

The water tunnel test section has a square cross-section. However, as a simpli- 
fication, the grid was made axisymmetric by revolving the grid around the X-axis. 
The grid outer diameter was selected so that the tunnels walls were represented with 
a cylinder of the same cross-sectional area as that of the test section. This match- 
ing of cross-sectional areas captures the tunnel walls interference effect as described 
in [13]. Figure 7-7 shows one meridional plane of the RANS grid. To form the three- 
dimensional grid, the meridional plane is rotated about the X-axis. The chosen grid 
is of size 89 axially by 33 radially by 144 azimuthally. Due to computational limita- 
tions, the grid is not refined enough at the inner and outer boundaries to capture the 
boundary layer. 
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7.2.2 Rotor Experimental Comparisons 



The computed force coefficients are shown in figure 7-8. The forces are shown ver- 
sus number of iterations between the lifting-surface computation and the viscous 
flow solver. The thrust and torque computed in the coupled analysis using the 
PUF14/DTNS3D method are within 7.3% and 5.8% respectively, of those obtained 
experimentally. The same case, save with a finer RANS grid, has been checked using 
PBD/DTNSax with thrust and torque computed within 3% and 1.5% respectively, of 
those obtained experimentally. It is suspected that the much of the difference between 
PUF/DTNS3D and PBD/DTNSax results is due to the same effective inflow error 
discussed in section 7.1. No data on the stator is available when operating upstream 
of the rotor. 

Adjustment to the calculated rotor thrust due to the presence of a hub vortex was 



MIT Pre-swirl 
grid and blades 




-3 -2.5 -2 -1.5 -1 -0.5 
Upstream 



4.5 5 5.5 6 

Downstream 



Figure 7-7: One meridional plane of the RANS grid is shown with the stator and rotor blades 
overlayed. The three-dimensional RANS grid is formed by revolving this meridional plane about the 
X-axis. 
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Validation of PUF-14/DTNS3D using MIT Pre-swirl case 




Lifting-Surface/RANS Coupling Iteration Number 

Figure 7-8: Computed MIT Pre-swirl force coefficients versus coupling iterations. 

calculated by the method of Wang [40] as shown in [7]. The estimated hub vortex axial 
force was less than 0.2% of the calculated rotor thrust and was considered negligible. 
The change to the calculated rotor thrust due to hub gap forces, as described in [33], 
was likewise negligible. 

7.2.3 Stator Secondary Flows 

The stator is represented in RANS as bod}' forces at the RANS cells that correspond 
to the stator blade positions in the lifting-surface code. The flow field responds to the 
RANS-stator body forces in a similar manner as if a physical blade was present in the 
flow. Figures 7-9 and 7-10 show the flow field at X/R rot0 r = —0.62, about mid-chord 
of the stator blade. The flow accelerates on the suction side of the blade and slows 
on the pressure side. Near the tip, the flow tends to roll over from the pressure side 
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to the suction side subsequently forming a tip vortex. 

The major goal of this thesis is to capture the one-per-revolution forces in order to 
predict maneuvering forces. The stator secondary-flow detail will contribute to the 
rotor blade-rate forces, not the once-per-revolution forces. While the stator secondary 
flows are of interest, they are not the center of effort for this thesis. Further studies 
of the stator flow details are left to future work. 

7.2.4 RANS-Stator Thickness and Non-dimensional Time 

The non-dimensional time variable Ng controls the number of stator force positions 
written for conversion into RANS body forces. The stator force position that corre- 
spond to stator blade positions will have bod)' forces present; the other stator force 
positions have forces which have been set to zero. 3 Figure 7-11 illustrates the stator 
force placement in RANS. 

The RANS-stator ‘'thickness” is defined as the azimuthal extent of those RANS 
cells where the forces from a blade overlap the RANS cells. Overlapped RANS cells 
approximate the stator shape. The approximation improves with smaller RANS cells 
and higher Ng. The thickness of the RANS-stator approximation of the stator blade 
will strongly affect the local flows around and behind the stator. A RANS-stator 
approximation with smaller RANS cells leads to more defined secondary flows sur- 
rounding the stator. 

When N stator Blades — Ng , each stator must occupy the entire volume associated 
with its sector so that the stator forces are distributed in the full circumferential 
volume. In a uniform inflow, the RANS body forces and the resulting flow field 
would appear circumferential mean. However, in the 3D coupling, the apparently 
circumferential-mean flow field would not be an accurate representation of the physical 
problem. 

3 The program, PUF-I4, is written so that a stator BVP is calculated just the same as a rotor - except of coarse, 
that the stator does not rotate. As little stator-specific logic as possible was used in PUF-14. Doing this may make 
the computation time longer than necessary for a stator, but the payoff is a less complicated methodology. Since the 
stator is a special case of a rotor, i.e. not rotating, the variable Nq does take on a second meaning. 
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Recall that the forces transmitted to RANS are the time-average forces. In a broad 
sense, these time-average forces are generally obtained as the product of the solved 
blade circulation and the RANS calculated time-averaged total velocity. Normally, 
a RANS total velocity is used to obtain the effective inflow as the difference of total 
velocity minus time-average induced velocity. However, for a stator, the time-average 
induced velocity is equal to the local induced velocity. Thus, the RANS body forces, 
which would be distributed circumferentially when NstatorBiades = Ng, are not consis- 
tent with the local induced velocity. The corresponding lifting-surface/RANS solution 
would be in error. 

The error due to inconsistency is removed by increasing Ng so that the angular 
width of the stator body forces is small. Small angular width allows the formulation 
to be more consistent by causing local effects to become more pronounced. Based 
on the MIT Pre-swirl case (which has nine stator blades), a 10° width distribution 
of body forces in RANS seems reasonable. A 10° RANS-stator width corresponds 
to N$ = 72 4 and is seen in figure 7-12 as the abscissa of 0.1. Interpolations of the 
rotor mean harmonic force convergence plot indicates 1% error in convergence occurs 
at about an Ng = 60. Consequently, a minimum Ng of 60 should be used for this 
stator. Obviously, convergence of higher harmonic forces requires thinner RANS- 
stator width, which requires Ng to be higher. For instance, for this case, a Ng of 60 
would yield about a 1% error in the mean force and about a 10% error in the 9 th 
harmonic force. 5 Additionally, for cases with higher numbers of stator blades, it is 
expected that thinner RANS-stator width is required to maintain proper local effects. 
With experience in stator validation, perhaps the proper setting for the RANS-stator 

4 360°/72 = 5°, but the 1 inear interpolation algorithm places the zero forces at the adjacent angular locations so 
that the RANS-stator width becomes 2x5° = 10°. 

5 This methodology may have difficulty if Ng is very large and the RANS cells are very small in the tangential 
direction. The difficulty arises because the RANS-stator width becomes very small (the stator forces occupy an 
unrealistically thin volume of RANS cells). As a result, the stator forces are introduced into too small a volume in the 
RANS domain, which leads to high gradients in the induced velocities. The high gradients would lead to artificially 
high dissipation near and downstream of the stator. This high dissipation may wash much of the stator force effect 
out of the flow field prematurely. In the case of a stator upstream of a rotor, the wake dents from the upstream stator 
may get too thin to propagate properly to the downstream rotor. As a result, the higher-harmonic rotor forces due 
to the stator’s effect on the velocity field may begin to be incorrect. 
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width can be based on the ratio of RANS-stator width to number of stator blades, or 
somehow related to an estimated momentum thickness. 

Physically, the stator induces very local effects such as those shown in figures 7- 
9 and 7-10. Certainly, it is possible to represent the global influence of the stator 
reasonably well using circumferential-mean lifting-surface/RANS coupling. To switch 
to circumferential-mean stator representation is a fairly straightforward adjustment 
to the methodology. 6 However, the present methodology represents the stator using 
local influences to capture the once-per-revolution variations in stator loading which 
will be present during a maneuver. 



6 To represent the stator as a circumferential-mean blade row, two areas of change are required. First, the body 
forces must be smeared circumferentially. For instance, by adjusting the variable N e to equal NstatorBladcs • Second, 
the time-averaged induced velocity must be calculated using circumferential-averaging. 
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Figure 7-9: Local axial velocity formed by the stator body forces. 




Figure 7-10: In-plane velocity formed by the stator body forces. 
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Figure 7-11: RANS-stator width illustration 
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Figure 7-12: Convergence of Kt and I\q for mean and 9th harmonic rotor forces verses RANS-stator 
width. 
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Chapter 8 

Coupled Analysis of an 
Open-Propeller Slender Body 

8.1 Huang Body 1 Description 

This test case is one of a series of axisymmetric bodies which was tested by Huang with 
and without a propeller present [IS]. The same forebody (DTMB Model 5225) was 
used for all experiments while a series of afterbodies with increasing conicity were 
examined. Boundary layer profiles, skin friction, pressure tap and propeller force 
measurements were obtained and used by Huang to perform research in the areas 
of propeller/hull interaction and turbulence modeling [16, 17, 38]. These afterbodies 
have been used extensively for validating solutions to the effective wake problem using 
analytic and numerical methods [8, 36, 44, 45]. 

The afterbody considered here, Afterbody 1, is a non-separating stern with a low 
tailcone angle. A profile view of the entire body (Model 5225-1) is shown in figure 
8-1. It was tested in the presence of an open rotor in wind tunnel and towing tank 
facilities. The Reynolds number based on body length was 5.9 million in the wind 
tunnel tests that were used for comparison here. 

Propeller model 4577 is a 6.0-inch diameter, seven bladed, wake-adapted, alu- 
minum propeller. The geometry is shown pictorially in figure 8-2. Open-water per- 
formance data for this propeller was obtained from propeller 4567A using a propeller 
boat in a towing tank. Propeller 4567A is a 11.9 inch diameter duplicate of propeller 
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Figure 8-1: Pictorial representation of Huang Body 1 (DTMB Model 5225-1). 



4577. 




Figure 8-2: Pictorial representation of Propeller 4567A. 



8.2 Lifting-Surface/RANS Coupling Example 

8.2.1 Nominal Flow 

The nominal flow was calculated using a three-block RANS grid. The number of cells 
axially, radially and azimuthally was: 97x57x33, 73x57x33 and 73x57x33. The axial 
number, radial number and their distribution were selected based on convergence of 
the nominal inflow. The number of cells azimuthally was selected to provide the 
largest three-dimensional grid size that could be run on the computers at MIT. The 
results for all yaw angles except those near 0° are, most likely, not well converged with 
respect to the RANS grid; however, the trends around 0° remain valid and illustrate 
the advantages of three-dimensional coupling using PUF-14. 

Figure 8-3 shows the nominal inflow at xj L = 0.977. The serial-UNCLE results 
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appear to over-estimate the axial velocity over most of the propeller span. The 
axisymmetric boundary layer profiles were taken from Black [1]. 




UJV- 



Figure 8-3: Nominal inflow on Body 1 at x/L = 0.977 from experiment, axisymmetric DTNS and 
serial-UNCLE models. 



Obtaining the Grid 

The grid is obtained using inmesh [6] to form one meridional plane, which is composed 
of three zones. Then, the meridional plane is rotated about the X-axis to form the 3D 
grid. Figure 8-4 shows the overlap of the blade outline and the meridional plane. As 
a minimum , at least 10 RANS cells should evenly overlap the blade outline in both 
spanwise and chordwise directions. The figure shown here illustrates the absolute 
minimum RANS-cell/swept-blade overlay. 

Good overlap of the blade grid by the RANS grid is necessary for the following 
reason. The lifting-surface method calculates the time-average velocity at each control 
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point. The calculation is done using potential flow analysis and is therefore exact for 
the calculated blade loading. Large gradients in the time-average velocities could exist 
across the blade chord. To calculate the effective inflow, the time-average velocity 
is subtracted from the RANS total inflow. If the RANS flow field does not contain 
sufficient cells in way of the blade to capture the streamwise velocity gradients, then 
the results from the subtraction will be incorrect. 




Figure 8-4: Overlay of the blade swept outline and the RANS grid. 



8.3 Methodology of Coupling with RANS 

The lifting-surface methodology uses the total inflow extracted from RANS to calcu- 
late the effective inflow, then to solve the boundary value problem. Figure 8-5 show 
the converged total inflow for Huang body 1 at a 30° yaw to port, or 30° drift angle. 

To allow the entire RANS flow field to respond to the propeller’s influence, the 
time-average lifting-surface forces must be placed into the RANS domain. The time- 
average axial force is shown in figure 8-6 at every 3 rd position. The contour variable 
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Figure 8-5: This figure shows the total inflow at every 3 rd velocity disk. 



is h'r / ( '/J * Area). 

The time-average velocities provide the feedback to the PUF code so that the ef- 
fective inflow can be determined. By subtracting the previous time-average velocities 
from the current total velocities from RANS, PUF-14 calculates the effective inflow, 
which it uses to solve the boundary value problem at each timestep. Figure 8-7 shows 
the time-average axial velocity. 

8.4 Converged Coupled Results 



The results of the coupled lifting-surface/RANS method will produce the full flow 
field including the effects of the propeller on the flow field. Additionally, the PUF-14 
analysis yields the blade-row specifics such as blade pressure, local-blade forces and 
blade-row forces. Figure 8-8 shows the convergence of the mean I\'t and I\q at zero 
drift angle. 1 

*Each iteration took approximately eight hours on a DEC Alpha 333 MHz workstation. Of the eight hours, 
about 35 minutes were used to execute PUF-14 and the remainder was used for the execution of serial-UNCLE. The 
interpolation-weight setup program executed only every third coupling iteration and required about 40 minutes to 
execute. 
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Figure 8-C: Plot of the time-average axial force. 




Figure 8-7: Plot of the time-average axial velocity which is used to determine the effective inflow. 



Table 8.1 shows the experimental results compared to the calculated results. The 



first four rows are repeated from the work of Black [1]. These rows show the Propeller 



4577 performance in the experimental nominal inflow and that computed using differ- 



ent turbulence models when coupling PBD-14/DTNSax. The PUF-14/serial-UNCLE 
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results are shown in the last row. The serial-UNCLE code uses a Baldwin- Lomax al- 
gebraic turbulence model which gives a nominal inflow that resembles the DTNSax 
Baldwin-Lomax results at the hub and the k-e results from about mid-span outward. 
The results for the different turbulence models are compared in figure 8-3. Since 
the nominal inflow most closely matches the k-e DTNSax results over much of the 
blade span, it seems reasonable that the PUF-14 calculation of A'r and I\q should 
also match the k-e DTNSax results. In fact, since the serial-UNCLE seems to over 
predict the axial velocity compared to the k-e, it is expected that the propeller forces 
calculated from that inflow will be lower than expected from the k-e alone. Therefore 
the PUF-14/serial-UNCLE results are expected to be similar or slightly lower than 
the k-e DTNSax results for I\j. This is, in fact, the results that are obtained, and 
this supports the need for an improved turbulence model in serial-UNCLE. 2 





I\T 


error 


K q 


error 


Experimental inflow 


0.2902 


- 


0.05367 


- 


PBD-14/DTNSax k-e 


0.2520 


-13.2% 


0.04751 


-11.5% 


PBD-14/DTNSax Baldwin-Lomax 


0.2716 


-6.4% 


0.05081 


-5.3% 


PBD-14/DTNSax Modified Baldwin-Lomax 


0.2876 


-0.9% 


0.05327 


-0.7% 


PUF-14/serial-UNCLE Baldwin-Lomax 


0.2403 


-17.2% 


0.04751 


-11.5% 



Table 8.1: Propeller 4577 performance at zero drift angle using various turbulence models. 



To illustrate the prediction of maneuvering forces, the body was subjected to 
several drift angles. This simulates that the body is in a yaw, pitch or in the steady 
rotation during a maneuver. The mean A’r and A q are plotted against the drift angle 
in figure 8-9. Note that due to computational and time limitations, only 33 RANS 
cells were used in the azimuthal direction around the body. Thus, the RANS grid 
may not be fine enough to correctly capture the physical behavior of the flow field at 
the larger drift angles. 

Figure 8-10 shows the Huang body 1 at a 30° yaw to port. Streamlines show the 
general direction of the flow. The axial propeller force is shown in the contour plot 

2 Note that the RANS grid used in PUF-14/serial-UNCLE coupling is different from the grid used by Black. Slight 
differences could be due to grid issues. 
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on the blade surfaces. Evidence of the interaction of the propeller forces on the body 
RANS solution can be seen by the trajectory of the streamlines passing near the body 
stern. 

Figure 8-11 shows the trends in the harmonic amplitudes of the forces on a blade. 
The forces are shown in the xyz-coordinate system that is rotating with the blade. 
Similar trends exist for the x, y and z moments on a blade. Figure 8-12 shows the 
plot of XYZ-shaft-and-bearing forces versus rotation of the propeller in the inertial 
reference system fixed to the body. Similar trends exist for the X, Y and Z moments 
on the shaft. 

8.5 Maneuvering Forces 

The coupled lifting-surface/RANS method can be used to calculate the forces on a 
body during a maneuver. The body forces and moments are obtained by integrating 
the pressure and drag forces on the body from the RANS solution, and by adding 
the shaft-and-bearing forces from the PUF-14 solution. Since the RANS solution 
includes the presence of the propulsor, the RANS forces and moments reflect the 
actual thrust required to move the body at the specified speed (i.e. thrust deduction 
is included). The propeller forces are obtained from the time-domain lifting-surface 
solution. Figure 8-13 shows the coordinate systems of the lifting-surface analysis 
(PUF) and the RANS solution (RANS) relative to the body-centered coordinate 
system (MAN) used for maneuvering force descriptions. 

Figures 8-14 and 8-15 show maneuvering data of lateral forces and yawing moment, 
respectively. The “Nominal Experiment” is the measured data without the action 
of the propulsor. The experimental nominal data has a slight bias which can be 
seen by the non-zero value at zero degrees. This data can be calculated in RANS 
by converging the RANS nominal solution (i.e. not coupled with PUF-14). The 
calculated data from RANS alone is labeled “Nominal RANS”. The coupled results, 
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“Propelled PUF-14/RANS”, show the maneuvering forces including the influence of 
the propeller. No propelled experimental data was available for this case. Note that 
the RANS grid had only 32 cells azimuthally and will exhibit significant errors when 
the yaw angle becomes large. Calculated data greater than about 6 degrees becomes 
erroneous due to RANS grid azimuthal resolution. 

The propeller contributes to maneuvering forces in two ways. First, the propeller 
shaft-and-bearings experience forces due to the interaction of the blades with the 
incoming flow. The shaft-and-bearing forces directly place forces on the body. Second, 
the propeller interacts with the nearby flow field to change the shear and pressure 
forces on the body, i.e. thrust deduction. 

Without a propeller present, the body in yaw will experience no significant heave 
force nor pitch moment. However, with a propeller operating, both these components 
are present. Figures 8-16 and 8-17 show this force and moment, respectively. While 
the magnitudes of these are smaller than those shown in the preceding figures, these 
may have an important effect when acting over a long time during a maneuver. The 
secondary force and moment originate from the steady shaft force as seen in figure 8- 
12. 
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Figure 8-8: Convergence of AY and Kq versus lifting-surface/RANS iteration for Huang body 1 at 
zero drift angle. 




Figure 8-9: Converged AY and Kq versus drift angle for Huang body 1. 
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Figure 8-10: Huang body 1 coupled lifting-surface/RANS solution (Color reproduction in figure B-2.) 
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Figure 8-11: Force trends on a blade in the blade-coordinate system associated with varying the drift 
angle. 
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Figure 8-12: Force trends on the shaft in the inertial coordinate system associated with the varying 
drift angles. 
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Figure 8-13: Coordinate system depiction 
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Figure 8-14: Lateral force on Huang body 1 




Figure 8-15: Yawing moment on Huang body 1 
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Figure 8-16: Heave force when in yaw. 
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Figure 8-17: Pitch moment when in yaw. 



103 



Chapter 9 

Performance Prediction for a 
Multi-Stage, Ducted Propulsor 

9.1 Introduction 

The ducted pre-swirl unit initially designed by the Black [2, 28] was redesigned using 
current design tools during Engineers degree research by the author [41]. A description 
of the design procedure and exact geometry is presented by Warren [30. 41]. This 
geometry is to be used as a notional integrated stern that can be analyzed by other 
investigators. While experimental data for this geometry is not available, calculations 
for this geometry will be used to verify the method developed herein for an integrated 
stern. A pictorial representation of the geometry, named the 5/ren/an, is shown in 
figure 1-1. 

The diameter based Reynolds number for the stator and rotor were 2.3 and 1.9 
million, respectively. The RANS grid used here was derived from the one used by 
the author with modifications to reduce the computational time and memory require- 
ments of the serial-UNCLE solver. The current grid has good discretization axially 
and well as radially. The y+ is maintained between 1 to 4 for the first cell off the 
body. Azimuthally, the grid is relatively coarse, using only 33 cells to cover the full 
360° arc. As such, these coupled results are not expected to be accurate at large 
angles of attack. The azimuthal limit of 33 cells was necessary due to computational 
limitations on the computers used for this study. 
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It should also be noted that 33 cells were chosen to be an integer multiple of the 11 
stator blades. This choice reduces harmonic errors, which may have been introduced 
by the interpolation of the stator-body forces onto the RANS grid. Additionally, 
to help compensate for the coarse azimuthal discretization, the stator B-spline was 
rotated by 5° about the x-axis. 1 This rotation places each entire stator blade in one 
RANS azimuthal set of cells. Doing so will concentrate the body force and thereby 
strengthen the local velocity gradients. Higher velocity gradients should minimize the 
errors due to the different induced velocities calculated in the RANS domain and in the 
PUF domain. 2 Finally, the stator boundary value problem uses the circumferential- 
mean inflow. This is necessitated by time constraints. Future validation should use 
the full 3-D inflow to obtain stator side forces and moments. 

9.2 Example Case 

The case is the first verification of the whole methodology developed throughout this 
thesis bundled into one example. The case is run at drift angles of 0°,2 o ,4° and 6°. 
The total inflow is interpolated onto a set up non-radial conical surfaces. The cone 
surfaces provide velocity information in the region from just upstream of the blade- 
row to a few blade-row radii downstream. Figure 9-1 shows the cut-away volume 
representation of the cone surfaces for the stator in 4 ° drift angle. The stator blade 
outlines are shown in the figure. 

The body stern is shown in figure 9-2 with a cut-away view of the duct. The stator 
and rotor transition wakes are also shown. The wake shapes conform to the body 
geometry, subject to boundary layer interactions. The wakes do not follow cylindrical 
geometry assumptions that are present in previous unsteady lifting-surface methods. 
These examples of the wakes illustrate the benefit of wake-adapted modeling. Addi- 
tionally, the wake interaction does not take place using potential flow singularities; 

Alternately, the RANS grid could have been rotated. 

2 The stator B-spline rotation is not necessary, but should partly compensate for the coarse RANS azimuthal 
representation. If the azimuthal number of cells was much greater, say 110, then this rotation would not be necessary. 
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Figure 9-1: This figure shows the volume used to perform the stator analysis. This case is at 4 ° drift 
angle. Tangential velocity is shown in the contour (Color reproduction in figure B-3). 



instead the body force modeling in the viscous solution accounts for the interaction. 
The stator blade-row is analyzed using potential flow techniques in a specified inflow 
such as figure 9-1, then, separately, the rotor is analyzed in its specified inflow. Thus, 
each blade-row analysis is performed independently and singularities from different 
blade-rows do not explicitly interact. Both the stator wake and the rotor wake follow 
the circumferential-mean total flow field. On average the wakes should be force free; 
however, no explicit effort enforces them to be force free. Of course, the assumption 
of a force free wake also neglects secondary effects such as wake defects from blade 
boundary layers, cascading blockage effects from downstream blade-rows, etc . and is 
subject to the assumption of time-average viscous modeling. 

The serial-UNCLE forces and moments are calculated by integrating the pressure 
and shear stress on the body and duct. The stator and rotor forces and moments are 
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Stator and Wake 




Figure 9-2: Wake-adapted grids interact through the viscous solver. 

calculated by the lifting-surface method. Figure 9-3 shows the lateral force density 
contours on the body, duct and blade-rows. The body contours are the lateral force 
in the RANS coordinate system shown in the figure. The contours on the blades are 
in the blade-centered coordinate system, which rotates with the blade (be. contours 
are the tangential force). All contours are non-dimensionalized by the body radii to 
be consistent with the serial-UNCLE output. This figure illustrates the lateral force 
density for a body at a 4 ° drift angle from the starboard bow. The stagnation point 
on the body has shifted to starboard and, in general, the starboard side of the body 
feels a force directed to port. The stern shoulder has a strong force to starboard that 
is generated as the flow accelerates at the shoulder due to both the stern contraction 
and the yawed inflow. 

The flow that is a few propeller radii away from the propulsor is not strongly 
influence by the propulsor. However, in the propulsor region, the forces on the body 
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are influenced by the stator and rotor interactions. The stator body forces are at 
discrete locations in the RANS domain, admittedly with a somewhat broad angular 
representation. At the stator root, the triangular contour patches of higher lateral 
forces indicate the discrete stator interactions. That is, forces are present where the 
stator blade overlaps the RANS cells and no forces are present in the void in between 
the stator blade. The rotor forces, being time- averaged, fill the entire swept volume 
of the rotor blade-row with forces that smoothly vary in all spatial coordinates. The 
rotor body forces influence the RANS solution with an axial pressure rise. The rotor 
influence is most clearly seen by examining the inside of the duct. On the inner port 
side on the duct, the lateral force varies from strong negative to strong positive. This 
arises due to the pressure drop due to the stator and then a pressure rise due to the 
rotor. 
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Figure 9-3: Force contours on the body, duct and both blade-rows (Color reproduction in figure B-4.) 
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During the course of this work and other related work, one major concern was 
noted associated with the viscous modeling. For some cases, such as this Sirenian 
body but not the Huang Body 1 present in chapter 8, the RANS viscous solution 
oscillates around a mean value. The amplitude of oscillation does not decrease with 
more RANS iterations and the residuals do not improve. 

The converged-coupled solution for Sirenian at 4° drift angle was run until the ro- 
tor thrust and torque converged. The RANS solution had also converged as indicated 
by no further reduction in the residuals. However, the RANS solution had not actu- 
ally converged to a stable solution, instead the RANS solution was oscillating within 
a narrow range. More iterations of the RANS solution did not stabilize the solution; 
more iterations continued to oscillate the solution about some mean value. Presum- 
ably, the center of oscillation is the correct solution. The dotted lines in figure 9-4 
show the lateral force integrated around each body section from RANS solutions that 
were several hundred iterations apart. The bold line in an average of all solutions 
spanning between the dotted lines. The dotted lines show some typical extremes in 
the variations. 

In addition to the oscillation of the solution, it is surprising that the lateral force 
integration does not vary smoothly. For this Sirenian body, the geometry varies 
smoothly from the bow to the stern. The duct is the only discontinuity. Presumably, 
if many more force predictions (dotted lines) were averaged, then the average lateral 
force (solid line) would be smoother. It is troubling to find the solution by averaging 
many “solutions”. Interestingly, the force integration over the body stern does not 
seem to oscillate. The observed oscillation could simply be due to poor grid quality, 
or could indicate a more fundamental problem. Similar oscillations have been noted 
in axisymmetric RANS solvers on other body shapes. Admittedly, in the hands of a 
more accomplished RANS user, these oscillations may not be present. This concern 
is presented for completeness and is not meant to invalidate the viscous modeling in 
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RANS, but to simply point out an area of concern. It is recognized and supported 
that viscous modeling in RANS offers an extremely powerful tool for computational 
modeling. However, since the effort of this thesis in not directed toward the RANS 
modeling, research in the RANS modeling is left to other universities and research 
centers. This thesis effort centers on the lifting-surface blade-row modeling and the 
interaction with the given RANS solution. 




Figure 9-4: Lateral force ( Fz ) verses body position 



9.3 Maneuvering Forces 

Figures 9-5 and 9-6 show calculated maneuvering data of lateral forces and yawing 
moment, respectively. These forces and moments are non-dimensionalized by the 
body length. No experimental data are available. The trends follow those observed 
in the Huang Body 1 validation in chapter 8. The calculated data from RANS alone 
is labeled “Nominal RANS”. The nominal case is expected to be on the verge of 
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separating at the body stern due to the relatively full stern. While no evidence of 
separation was observed in the nominal solution, it remains important that the viscous 
solver be able to capture separation. At larger yaw angles, separation may be more 
likely. The calculated-coupled results, labeled “Propelled PUF-14/RANS” , show the 
maneuvering forces including the influence of the propeller. Note that the RANS grid 
used in this analysis had only 33 cells azimuthally and may exhibit significant errors 
the at larger yaw angles. 

The propulsor contributes to maneuvering forces in several ways. First, the rotor 
shaft-and-bearings experience forces due to the interaction of the blades with the 
incoming flow. The shaft-and-bearing forces directly place forces on the body. Second, 
the stator experiences forces that act directly on the body. Recall that, due to time 
constraints, this stator analysis is performed in the circumferential-mean inflow, which 
results in zero side forces. Still another contribution from the propulsor is a change 
to the nearby flow field which changes the shear and pressure forces on the body, i.e. 
thrust deduction. Final, the duct interacts with each of these modify both the forces 
on the components and on itself. The coupled lifting-surface/RANS method provides 
the mechanism to quantify their interaction. In this example, the propulsor effects 
on the maneuvering forces are relatively small, approximately 7-10% change from the 
nominal solution. 

Figures 9-7 and 9-8 show the heave force and pitch moment, respectively. The 
secondary force and moment originate from the steady shaft force and moment. In 
both these figures, some RANS error can be seen by observing that the nominal 
calculated results should have zero heave force and pitch moment. Like the Huang 
Body 1 case, the magnitudes of the heave force and pitch moment are smaller than 
lateral force and yaw moment, respectively. While the heave force and pitch moment 
are smaller, these may have an important effect when acting over a long time during a 
maneuver. Additionally, the propulsor-induced force and moment strongly depend on 
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the body-boundary layer interaction with the propulsor and may be relatively small 
for this example. 
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Figure 9-5: Lateral Force 




Figure 9-6: Yawing Moment 



113 



0.0003 


— A — Nominal RANS 
— hi Propelled PUF14/RANS 


0.0002 




0.0001 




£ o 


- 


- 0.0001 


• 


- 0.0002 




- 0.0003 


* » i ' ' * ' 1 ' 1 * LU 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ! 1 1 1 1 1 J 1 1 F * 

-7 -6 -5 -4 -3 - 2-10 1 

BETA 

Figure 9-7: Heave force when in yaw. 
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Figure 9-8: Pitch moment when in yaw. 
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Chapter 10 

Conclusion and Recommendations 



An original computer code, PUF-14, has been written to support the new methodol- 
ogy developed throughout this thesis. The lifting-surface method. PUF-14, is coupled 
with a RANS code provided by our ONR sponsor. The RANS code is used to obtain 
the inflow velocity field that the lifting-surface code uses to calculate the unsteady 
forces. Unsteady forces are generated due the rotating propeller in a spatially-varying 
inflow. Time-averaged, but spatially-varying body forces are introduced into a three- 
dimensional volume to represent propulsor stages in the RANS flow field. The entire 
RANS flow field responds to the blade-row presence. In turn, the RANS flow field is 
used again for the lifting-surface analysis of the blade-row. By alternately updating 
the lifting-surface and the RANS solutions, the blade-row forces and RANS flow field 
converge to the appropriate solution. 

The new treatment of unsteady force calculations should greatly improve propul- 
sor prediction capabilities. The new treatment is believed to be practical in both 
computational load and in representing the physical hydrodynamic characteristics of 
today’s complex propulsors. 

10.1 Improvements on the Current Method 

In developing this methodology, many possible improvements have become obvious 
and should be implemented given time and resources. Some improvements may make 
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only small improvements in accuracy; others may be more significant. The improve- 
ments are detailed below in no particular order. 

• Wake and blade lattice model which includes rotational variation 

The current method assumes the streamtubes which convect the shed and trailing 
vorticity are the circumferential-mean streamtubes. This assumption allows one 
set of influence functions to be used at all blade positions. An improvement would 
model the rotational variations in the wake lattice so that the wake follows the 
full three-dimensional velocities. 

• Stand-alone unsteady wake alignment 

Currently, to align the wake when not coupling with RANS, convection velocities 
are used to stretch the wake to account for induced velocity. Therefore, it may 
not be free of forces as it should be. A wake alignment routine would enhance 
the use of PUF-14 when used as a stand-alone analysis tool. 

• Trailing edge sheet model 

More research is needed to extend the 2-D Lagrange Kutta condition to 3-D. 
The 2-D model seems very promising and. once expanded into 3-D. may greatly 
improve the accuracy of higher harmonic forces. 

• Tip gap model 

The current methodology does not attempt to account for real-fluid effects at 
the blade tips. A semi-porous tip gap model could be incorporated into this 
method. The model could be based on the circumferential-mean tip-gap effects 
or could encompass a time- varying nature. 

• Viscous-load and thickness-load coupling 

The growth of viscous boundary layers on the propulsor blades result in changes 
to the inviscid pressure distribution and generally a loss in lift. The displacement 
thickness of the boundary layer effectively changes the camber, thickness and 
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angle of attack of each blade section. Simularily, the thickness of the blades 
influence the blade boundary layer. The effects become more important for 
advanced-blade sections. The time-varying nature of the viscous effects would 
need consideration. 

• Blockage 

Blade thickness results in a blockage to the flow which affects the flow distribution 
past the propulsor and the duct, if one is present. The creation of the viscous 
boundary layer on the blade surface creates additional velocities that will have 
an effect on the through-flow and performance of downstream stages. These 
factors need to be considered if more accurate propulsor performance is to be 
predicted. 

• Improved algorithm for time-average induced velocities 

By far, the most calculationally intensive portion of the lifting-surface method is 
the calculation of the time-average induced velocities. This calculation requires 
the effects at every control point during one blade-interval passage from every 
vortex element on every blade and every sheet. A more advanced algorithm 
would improve the accuracy and greatly improve the speed of this routine. 

• Coupling with RANS using a single disk of velocities 

The guiding cases used in developing this methodology were multi-component 
propulsors with highly contracting stern flows. While this method works fine 
for open-water propellers, the additional computational load is a burden. Us- 
ability would be greatly improved for open-water propellers if the method were 
reduced to coupling with RANS via a single disk of velocities upstream of the 
propeller plane instead of the full 3-D volume of velocities. Since RANS veloci- 
ties would not be available downstream of the propeller to convect the wake, a 
wake alignment technique would be required to ensure the wake remains force 
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free. 



• Improved algorithm for tracking 2-D streamlines 

The current algorithm for tracing the 2-D streamlines in the circumferential- 
mean flow is not as robust as it could be. The algorithm struggles in the viscous 
sublayers and leads to difficulties when automating the coupling process for 
certain geometries. 

• RANS turbulence modeling 

The serial-UNCLE used for the latter work implemented a Baldwin-Lomax tur- 
bulence model. This was proven to be inadequate in the Huang body 1 compar- 
ison. More research should be directed to improve the turbulence models for the 
axisymmetric boundary layer flow and for the vortex-core flows associated with 
maneuvering bodies. 

• RANS parallel processing 

There is current research in RANS parallel processing algorithms. The method 
developed herein should be coupled to a parallel version of a RANS code. By 
far, the component in the coupled lifting-surface/RANS method requiring the 
most computer time and computer memory is the RANS code. Effort should 
continue to be directed at improving this component. 

10.2 Possibilities 

The stand-alone lifting-surface method may be useful in other applications. With 
some modifications, one possible use could be to evaluate the transitory response 
to step inflow changes. Both the transitory blade forces and the shaft-and-bearing 
transitory forces could be determined. Transitory response may be useful in designing 
the motor controller for a propeller. 

The coupled lifting-surface/RANS method could be used to model the appendages 
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of a body. For example, a control surface is not too different from a stator blade. Thus, 
the appendages could be modeled as “stator” blades; A four bladed “stator” would 
model four stern control surfaces. If the four blades were not oriented identically 
about the body centerline, then perhaps they could be modeled as four one-bladed 
stators. While the flow details would not be precisely correct, for a body at an angle 
of attack, the gross flow would be correctly captured. The coupled method would 
provide a relatively quick estimate of the maneuvering forces on an appended body, 
without ever having to remake a RANS grid. 

Finally, another possibility might use the coupled method interactively with a time- 
domain maneuvering simulation. The propulsor forces and moments could be calcu- 
lated at discrete time steps in a maneuvering simulation. The resulting body-forces 
could update the RANS solution, while the propulsor forces influence the trajectory 
of the body in the RANS domain. 

10.3 In Retrospect 

First, it must be acknowledged that neither PUF-14, the coupling methodology nor 
RANS codes are perfect. They all have shortcomings and all require careful attention 
to avoid the “garbage in - garbage out” syndrome. With effort, these methods can 
improve. 

The methodologies developed and incorporated into the stand-alone PUF-14 pro- 
vide the modern propulsor designer a tool to analyze trends in propulsor perfor- 
mance. Without RANS, the stand-alone methodology suffers by not possessing an 
automated wake alignment method. User input convection velocities fill this short- 
coming. However, as a propulsor analysis method, the less user input to the physical 
hydrodynamics the better the method. 

The coupled methodology consists of PUF-14 and a RANS solver. This method 
has great potential and should be invaluable to the modern propulsor designer. In 
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addition to trend analysis, the coupled methodology should provide relatively good 
agreement with experiment. The single-most important advantage is the ability to 
discretely model the component stages while avoiding the complex and computa- 
tionally intensive modeling of the rotating blades in the RANS domain. It is my 
sincerest hope that the methodology developed herein becomes invaluable to today’s 
and tomorrow’s propulsor designers. 
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Appendix A 



A Numerical Method for the 
Computation of Unsteady 
Propeller Forces 

Introduction 

During the computation of unsteady forces, many possible approaches could be adopted 
in locating the shed vortices within their corresponding time interval. Nearly any 
placement will converge to the correct answer given sufficiently small time steps. Ob- 
viously, the best numerical scheme is one that converges within the desired accuracy 
with minimal computations, i.e. the largest possible time step. 

In general, the spacing in the wake must be on the same order as the last spacing 
on the foil to get reasonable results near the trailing edge. Fine spacing near the 
trailing edge leads to the necessity of similarly fine spacing in the wake to remove 
the singular behavior at the trailing edge. Such fine spacing in the wake leads to 
considerable computational effort. 

Past unsteady vortex-lattice computational methods have used constant spacing 
on the foil with similar-sized spacing extended into the wake. Spacing arranged in 
this manner has provided good accuracy of the global forces with reasonable sized 
elements in the wake. 

How'ever, it is desired to resolve the gradients better on the foil, especially near 
the leading and trailing edges. A natural solution is to use Glauert constant angle 
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Figure A- 1 : Discrete vortices and control points along the chord. 

spacing, alias cosine spacing, on the foil. This choice places a fine grid of points near 
the leading and trailing edges, and a relatively coarse grid in the mid-chord region. 
Figure A-l compares constant and cosine spacing on a two-dimensional foil. 

Trailing Edge Singularity with Cosine Spacing on the Foil 

As shown in figure A-2, the flat-plate solution in a gust is sensitive to the placement 
of the vortex with the time-step element. Numerical studies such as Frydenlund and 
Kerwin [12] have examined the sheet strength as it transitions into the wake. The 
sheet strength should be continuous and smooth, except for a slope change at the 
trailing edge. 

The location of the discrete vortex within the first time step element in the wake 
causes dramatic changes in the sheet strength on the foil and the wake. A closer 
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Figure A-2: Trailing edge behavior of the vortex sheet strength. 



examination indicates that using an implicit Kutta condition allows the last control 
point to be in close proximity to the vortex in the wake. The closeness was avoided 
in previous formulations by using constant spacing on the foil and an explicit Kutta 
condition. 

Figure A-2 illustrates the strong dependence of the sheet strengths on the posi- 
tion of the first wake vortex and concisely illustrates the need for this formulation. 
The figure describes a linearized flat-plate in a sinusoidal gust at an arbitrary time 
during the cycle. The case with constant spacing on the foil serves as a guide to the 
correct behavior. All cosine-spaced cases have the identical spacing on the foil and 
wake except that the vortex placement within the wake time step element is varied. 
Although far from correct, placing the vortex at the quarter chord location within 
the element seems to provide the most reasonable compromise. Perhaps not quite so 
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coincidentally, the quarter chord placement of the vortex along with the three-quarter 
chord placement of the control point yields the exact solution in steady flow in the 
limit- of an infinite number of vortices. [20] 

This appendix discusses the formulation and implementation of a method to over- 
come the singular behavior of a discrete wake vortex near the foil control point. The 
formulation drastically reduces the computational effort while obtaining very similar 
accuracy. 

Formulation of the Two-Dimensional Problem 



Consider a two-dimensional thin hydrofoil advancing with constant speed U , which 
may be passing through a spatially varying velocity field. The linearized problem will 
be solved. A flat plat, with chord length, C, is situated on, or close to the interval 
(0, C) of the x-axis in a flow with speed U oriented in the positive x direction. 

The integral equation for the distribution of vorticity 7(.r ,£) over the foil may be 
written in the following form. 



,(,.«)+ / c 2«4 iC+ 

Jo x — L Jc x — n 



(A.l) 



£ Jc x — rj 

where £ and rj are dummy coordinates on the foil and wake, respectively, and ~, w is 
the strength of the shed vorticity in the wake. 

According, the Kelvin’s theorem, the total circulation of the system must remain 
zero, 

rC rUt 

/ 7(C»0 rf C+/ lw(ri,t)dri = 0 (A.2) 

Jo Jc 



so that the change in circulation on the foil must be followed by an equal change in 
circulation in the wake. 

The formulation is completed with a statement of the Kutta condition, which re- 
quires that 7(2, t) be bounded at the trailing edge. The Kutta condition can be made 
explicit by fixing the last bound vortex strength to be a value which satisfied the 
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desired behavior. Conversely, the Kutta condition can be implicit which is accom- 
plished by relative placement of the last control point on or very near to the trailing 
edge. 

Discrete Time Step Solution of the Boundary Value Problem 

For the discrete vortex model employed here, the governing integral equation, eq. A.l, 
is reduced to a system of linear algebraic equations: 

+ £(r.wn* = -K (A.3) 

j= 1 Ar=l 

The quantities are as follows: 

• N is the number of bound vortices 

• N w is the number of vortices in the wake 

• Bij are the influence functions which describe the induced velocity on control 
point i due to a unit strength vortex j located on the foil 

• Wi'k are the influence functions which describe the induced velocity on control 
point i due to a unit strength vortex k located in the wake 

• are the velocities acting at control point i due to the boundary conditions 

• (r&)jare the unknown bound vortices at position j on the foil at the current time 
step 

• (T s) k are the shed vortices at position k in the wake. 

The notation, r u s m , represents the discrete vortex strength in the time step el- 
ement. The superscript u in r u s m indicates that the circulation is counted in the 
un-subdivided intervals, which are counted with the index m. The numerical objec- 
tive of this model is to represent the r u ,$o with an alternate set of shed vortices, Fs k 
where k = 1 , 2,-*-,A r /. The vortices T u s m for m > 0 are left untouched and equal 
Fs k for k > Nj. 

The total circulation shed into the most recent time step element m = 0 in the 
wake is designated as T u s 0 . The T n represents the sum of the bound vortex strengths 
at time step n. Thus, Kelvin’s theorem can be discretized in equation A. 4 by writing: 

r u s 0 = —T n + T n _! (A. 4) 
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Profile View of a 2-D Foil 
and its wake 



Nntinnal Poil Outline 


• 


Foil Control Points 


r 


Foil Bound Vorticity 


o 


Wake Shed Vorticity 
using sub-divisions 


□ 


Wake Shed Vorticity 
using no sub-divisions 




Figure A-3: Pictorial of subdivisions in the first time element 

T, = |i;(r 4 )i I (A. 5) 

This new formulation of the linearized two-dimensional foil differs from previous for- 
mulations in that the most recent time step element is subdivided into Nj intervals. 
Each interval is represented by a discrete vortex. Thus, the shed vortices are sepa- 
rated into those within the most recent time step, which are unknown quantities, and 
into those from prior time step solutions which are known. Rewriting equation A. 3, 
we get: 

+ E(r«)tiVijt = -v { - £ (a.6) 

j=l k—\ k=N f 1 

where (Fs)/,. are known for k > N j but unknown for k <= Nj. 



132 



In order to maintain finite loading at the trailing edge of the foil, the vortex sheet 



must be continuous from the foil into the wake. Thus, by the relationship 



■<—hw l w - x u )) 



(A.7) 



we would like for the discrete vortex to be at least of order two so the sheet will 
be continuous. One simple way is to assume the T(^) is a polynomial of order A 
and use Lagrange interpolation to represent the discrete strengths in the wake. The 
interpolation essentially defines the distribution of vorticity in the subdivided interval. 
The total vorticity in the subdivided interval still remains as the only unknown in the 
wake. Since the total vorticity was unknown in the original 2-D formulation, using 
Lagrange interpolation does not introduce any more unknowns to the boundary value 
problem. Instead, the interpolation modifies the influence functions on both sides of 
the equation. 

The sheet vortex strength, multiplying the uniform shed wake element size U5t, 
can now be represented by Lagrange polynomials as follows. If r is the fractional 
time, or fractional distance, back from the present then 



A 



r Ain) = Y, L rn(T k )rs m 

m=0 


(A.S) 


Lm(r k ) = fl 


(A. 9) 



where represents the “product of”, A is the order of the Lagrange interpolation and 
m is the index which counts the un-subdivided intervals from the current “zeroth” 
time step. 

In the most recent time step element (in = 0), the discrete subdivision strengths 
can be obtained by applying Lagrange interpolation. Care must be maintained to 
ensure the total discrete strengths correctly model the integrated sheet strength. 
Equation A.S and A. 10 can be used to describe the unknown subdivision strengths as 
a combination of the current unknown strength and the known strengths of the prior 
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time steps. 



rs k = r A (r,) 



USt/N, 

USt 



(A. 10) 



For example, for a second order polynomial the Lagrange interpolation coefficients 
and the resulting equation describing the subdivision strengths are: 



r s k = (pkT u s 0 + qkF u Si + r k r u s 2 )— ^ j- f - with k = 1, • • • , N f 



Pk = 



t\ - 3 T k + 2 



qu = ~rl + 2 T k 



Tk = 



~k ~ n 



( A.l 1) 

(A. 12) 
(A. 13) 
(A. 14) 



It should be clear that r u Si = + i and that r u s 2 = r.s/v /+2 , and so on. which 

are all known quantities of previous solutions. The remaining unknown is r u ^ 0 which 
is obtained from Kelvin's theorem written as equation A. 4. After substitution and 
collection of the unknown quantities of the left hand side of the equation, the formu- 
lation with subdivision can be written as: 



SOU 

j=i 



1 



- yvT E ( n 



n - Tj 



Wi 



J k= i \ 7=1 r ° T i , 



l i,k 



= -Vi 



(A. 15) 



1 



*7 



N- ^ 



/ k=i 



n 



Tk ~ t j 



;=i r o - v 



Tn-l +e n 



Tk ~ Tj 



Ys 



m=l \j = 0 j^m 



Nf +m 



N w 



Wi, k - J2 (r s) k \Vi,k 

k=N, + l 



Or, for a second order polynomial, equation A. 15 simplifies to become: 



jv 

J = 1 

N 



N f 



ft. - TT E pMj, 



A 7 i=i 



= -v; 



(A. 16) 



J /V/ J\J W 

— t — [pkT n ~ i + ^-rs^+i + r fc rs^ /+2 ] w i'k - ^ 



Nf k= 1 " *=7^ + 1 

where p, and r are defined in equations A. 12 through A. 14. 

Two-Dimensional Results with a Second Order Interpolation 

Like figure A-2, figure A-4 shows a arbitrary time step for two cases. The first case 
is the more typical single vortex representation for each time step element. The 
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Figure A-4: Trailing edge behavior of the vortex sheet strength. 

discontinuity exists at the trailing edge due to the uneven spacing between the cosine 
spaced foil elements and the constant spaced wake elements. The second case shown 
in figure A-4 is using the newly developed formulation discussed earlier. The new 
formulation solves the boundary value problem by representing the first wake element 
more like a sheet as opposed to a single concentrated vortex. Sheet vorticity further 
away from the foil a sufficiently represented as a single concentrated vortex. However, 
near the trailing edge, a single vortex representation of the time step element causes 
erroneous results due to the singular nature of a discrete vortex approaching the 
trailing edge of the foil. 

The major advantage to this formulation is that it allows relatively coarse time 
steps while getting the accuracy obtained with much finer time steps. Figures A-5 
and A-6 illustrate this advantage by examining the foil vortex sheet strength at a 
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particular x/C location as the gust advances over the foil. In these figures, the first 
number in the legend is the time step element size relative to one another. The sec- 
ond number in the legend is the number of elements in the first time step (A/). For 
example, figure A-6 shows a case, labeled “2*, 1 element”, with no subdivisions and a 
very small time step. The relatively small time step is required to get the solution to 
converge to the correct results near the trailing edge. The corresponding subdivision 
case to comparison against is the ‘TOO*, 50 element”. With subdivisions, the case 
converges to nearly the identical results. This particular subdivision case uses 50 
subdivisions which divide the “100*” step into subdivisions of size “2*”. Thus, this 
case closely approximates the “2*, 1 element” case. The advantage is that the extra 
50 subdivisions are solved in the boundary value problem and are not required to be 
maintained in the wake past the first time step element. Additionally, the extra sub- 
divisions do not increase the number of unknowns. Instead, the subdivisions simply 
modify the coefficients of the simultaneous equations as shown in equation A. 15. The 
new formulation dramatically reduces computational effort while maintaining similar 
accuracy. 

Conclusion 

The new formulation dramatically reduces the computation effort to attain the same 
accuracy. The case studied herein applies a second-order interpolation; higher order 
interpolations should be evaluated to determine if they yield even more computation 
gain by increasing the time step size further. Additionally, this method should be 
applied to a three-dimensional unsteady code. The ultimate use for this method is 
to improve real life calculations in a three-dimensional unsteady method. 
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Wake Spacing Comparisions 
by observing the foil sheet strength 




Figure A-5: Vortex sheet strength behavior at x/C = 0.9 on the foil in a sinusoidal gust. Using a 
non-dimensional time step of 100* with Nj = 10 vortices in the first interval yields similar results 
as a time step of 10* and no subdivision. 
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Wake Spacing Comparisions 




Figure A -6: Vortex sheet strength behavior at xjC = 0.98 on the foil in a sinusoidal gust. The new 
formulation converges to similar accuracy while using significantly larger time elements in the wake. 
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Appendix B 

Color Figures 



PUF-14 Swept Volume Axial Forces 





m 



- 1.00 - 0.66 - 0.33 0.01 0.34 0.68 1.01 




Figure B-I: This figure shows a notional inflow with the corresponding time-average forces in the 
swept volume of the rotor (see section 6.3.3). 



139 









Underwater 
notional body 
and propeller 
at 30° yaw to 
port 

Color contours 
represent the 
computed 
axial force 



Figure B-2: Huang body 1 coupled lifting-surface/RANS solution (see section 8.4) 
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Figure B-3: This figure shows the volume used to perform the stator analysis. This case is at 4 ° 
drift angle. Tangential velocity is shown in the contour (see section 9.2) 
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Lateral Force ( F z ) 




- 0.2245 - 0.1677 - 0.1108 - 0.0540 0.0028 0.0596 0.1165 0.1733 



Figure B-4: Force contours on the body, duct and both blade-rows (see section 9.2) 
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